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STATIONARY  DISTRIBUTION  OF  PARTIAL  REMAINDERS  IN  S-R-T  DIGITAL  DIVISION 

Richard  Robert  Shively 
Department  of  Electrical  Engineering 
University  of  Illinois,  19&3 

The  execution  time  required  for  digital  division  in  electronic 
computers  can  be  reduced  if  a)  more  than  a  minimum  number  of  divisor  multiples 
are  available,  and  b)  provisions  exist  to  either  by-pass  or  shorten  steps  which 
generate  zero  quotient  digits.   A  division  algorithm  which  uses  this  redundancy 
to  increase  the  likelihood  of  zero  digits  (i.e.,  when  the  quotient  is  expressed 
in  the  recoded  form  implied  by  the  divisor  selections)  then  is  desirable. 

To  provide  a  basis  for  comparing  various  division  algorithms  in  this 
respect,  C.  V.  Freiman  observed,  using  a  Markov  chain  model  for  the  division 
process,  that  for  any  given  divisor,  D,  and  a  randomly  distributed  set  of 
dividends,  the  distribution  of  partial  remainders  at  successive  steps  approaches 
a  steady-state  or  stationary  distribution,  which  is  independent  of  the  distribu- 
tion of  dividends.   In  attempting  to  determine  this  stationary  distribution  as 
a  function  of  D  for  an  algorithm  known  as  S-R-T  division,  Freiman  found  that 
the  pointwise  solutions  provided  by  a  Markov  chain  are  inadequate  to  analyze  the 
function  for  certain  intervals  of  D. 

The  purposes  in  further  pursuing  this  work  started  by  Freiman  are: 

a)  to  obtain  a  more  complete  understanding  of  S-R-T  division, 

b )  to  derive  methods  which  might  be  found  applicable  in 
analyzing  other  division  algorithms,  and 

c)  to  provide  information  which  could  prove  useful  in 
formulating  new  division  algorithms. 

If  partial  remainders,  R.,  are  defined  to  be  in  one-to-one  corres- 
pondence with  quotient  digits  (rather  than  considering  only  those  associated 
with  nonzero  quotient  digits),  the  stationary  probability  density,  F(x),  which 


describes  R.,  can  be  shown  to  be  symmetric  with  respect  to  the  point  D,  F(D), 
for  all  D  in  the  allowed  range,  —  <  D  <  1.   This  symmetry  can  be  used  to  derive 
a  straightforward  method  for  evaluating  the  stationary  density  for  intervals  in 
the  divisor  range,  rather  than  for  individual  points „   The  method  essentially 
consists  of  locating  points  at  which  F(x)  is. discontinuous,  as  an  ordered 
process,  starting  with  points  of  discontinuity  in  the  relation  between  R.  and 
R.   ,  and  terminating  when  the  first  such  discontinuity  point  falls  in  the 
interval  (— ,  2D  -  p).   Then,  by  observing  that  the  jump  sizes  at  successive 
points  in  this  ordered  set  of  argument  values  form  a  geometric  progression,  and 
that  the  sum  of  these  jumps  is  — ,  the  solution  follows  directly.   Solution 
intervals,  i.e.,  neighborhoods  in  the  range  of  D  for  which  F(x)  has  similar 
form,  may  be  grouped  into  infinite  families,  based  on  the  sequence  of  operations 
used  to  locate  discontinuity  points,  and  when  considered  in  the  order  of 
increasing  complexity,  these  families  form  a  nested  pattern,, 

The  stationary  densities  describing  R.  can  be  used  to  determine  the 
shift  average,  <S>,  i.e.,  the  average  number  of  quotient  digits  formed  per  use 
of  the  adder.   <S>  is  equivalent  to  the  mean  recurrence  time  of  the  event 

"R.  >  — "  in  S-R-T  division,  and  is  therefore  the  reciprocal  of  the  probability 

1                                                   2D 
that  R.  >  p.   Using  properties  of  F(x),  <S>  can  be  shown  to  be  — =■  for 

0.75  <  D  <  1,  reaching  a  maximum  of  3.?  over  the  interval  0.6  <  D  <  0.75  •   In 

the  interval  0.5  <  D  <  0.6,  <S>  requires  an  infinite  number  of  sub intervals  to 

express  completely,  and  exhibits  a  nonmonotonic  behavior  at  certain  points. 

Analysis  of  other  division  algorithms  can  use  the  fact  that  a)  the 

piecewise  constant  functions  which  describe  the  stationary  density  can  be 

formed  by  locating  discontinuity  points  as  an  ordered  process,  and  b)  the  sizes 

of  jumps  at  successive  arguments  in  this  process  form  a  geometric  progression. 


1 .   BACKGROUND 

1.1   Introduction 

Numerous  factors  are  effecting  a  willingness  to  invest  in  redundant 
equipment  to  improve  digital  computer  operation  times.   Reliability  and 
decreasing  cost  of  components,  modularization  of  circuits,  and  miniaturization 
are  but  a  few  of  these.   In  arithmetic  processes,  some  commonly  considered 
applications  of  redundancy  within  the  structure  of  a  general  purpose  computer 
are:   a)  redundant  number  representation  to  reduce  or  eliminate  adder  carry 
propagation  time,  b)  multiple  shift  paths  to  reduce  the  number  of  steps  required 
in  iterative  processes  such  as  multiplication  and  division,  and  c)  generation 
of  more  than  the  minimum  number  of  multiples  of  the  operand,  again  to  improve 
multiplication  and  division. 

In  digital  division,  redundancy  in  the  generation  of  divisor  multiples 
may  be  used  to  serve  either  or  both  of  two  general  aims:   a)  providing  overlap 
in  the  partial  remainder  ranges  to  which  the  respective  divisor  multiples  are 
assigned,  and  b)  minimizing  the  expectation  value  of  the  partial  remainder 
magnitude.   The  first  of  these  is  particularly  important  when  signed  digits  or 
stored  carries  are  used,  since  distinguishing  between  nonoverlapping  intervals 
by  inspection  of  leading  digits  of  unassimilated  numbers  would  be  impossible. 
Reduction  in  the  partial  remainder  magnitude  can,  under  certain  circumstances, 
reduce  the  time  required  for  the  division  process.   (This  will  be  discussed  in 
Section  1.2. ) 

The  effectiveness  of  a  particular  digital  division  algorithm  (and  the 
redundant  equipment  it  may  require)  in  reducing  the  average  partial  remainder, 
R.,   can  be  measured  if  the  probability  distribution  (or  density)  functions 


Hereafter,  R.  will  be  used  to  denote  the  magnitude  of  the  partial  remainder 
occurring  on  the  ith  step  of  a  division  process. 
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which  describe  R.  are  known.   In  Section  1.3  it  will  be  shown,  using  the  theory 
of  finite  Markov  chains,  that  for  a  particular  divisor,  D,  the  distribution 
describing  R.  converges  with  increasing  i  toward  a  unique  steady-state  dis- 
tribution.  Characterizing  digital  division  as  a  discrete  Markov  process  is  due 
to  C.  V.  Freiman  (see  Ref.  4).   As  Freiman  discovered,  however,  the 
complicated  behavior  of  the  steady-state  distribution  of  R.,  for  certain  ranges 
of  the  divisor,  is  not  amenable  to  investigation  using  Markov  chain  techniques 
alone „ 

The  purpose  of  the  dissertation  is  to  derive  a  method  for  more  readily 
finding  the  steady-state  distribution  of  R.,  assuming  the  S-R-T  division 
algorithm  is  being  used,  and  to  apply  this  method  in  examining  the  function's 
behavior  in  those  ranges  of  D  where  detailed  analysis  was  previously  impractic- 
able.  Other  than  satisfying  a  curiosity,  the  worth  of  this  information  has  been 

2 
demonstrated  by  Professor  Gemot  Metze  who,  using  properties  indicated  in 

previous  efforts  and  proved  herein,  has  proposed  a  modification  to  S-R-T  divi- 
sion which  yields,  for  all  values  of  the  divisor,  D,  the  maximum  shift  average 
(i.e.,  maximum  number  of  quotient  digits  per  use  of  adder)  obtainable  under  the 
constraint  that  the  available  divisor  multiples  are  +D,  0,  and  -D. 

1.2   S-R-T  Division 

Descriptions  of  conventional  restoring  and  nonrestoring  division  can 

3 
be  found  in  various  texts  on  digital  computer  arithmetic. 

By  common  usage,  the  phrase  "S-R-T  division"  connotates  what  may  be 

regarded  as  nonrestoring,  binary  division  modified  to  utilize  the  facilities  of 


2 
'  See  Ref.  No.  7. 

3 

J  e.g.,  Ref.  No.  9. 
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a  floating-point  computer.   Sweeney,   Robertson,   and  Tocher   independently- 
observed  that  if  the  absolute  value  of  the  fractional  portion  of  the  divisor, 

7 
D,  is  restricted  to  the  normalized  range  .5  <  D  <  1,   then  any  steps  in  which 

the  partial  remainder  magnitude  is  found  to  have  leading  zeros  can  be  by-passed, 

and  zeros  inserted  as  the  corresponding  quotient  digits.   (After  being  recoded 

from  a  representation  using  +1,  0,  and  -1  as  digits  to  binary,  the  quotient 

digits  just  mentioned  will  be  all  zeros  or  all  ones,  depending  on  whether  the 

next  lower  weighted  nonzero  quotient  digit  is  positive  or  negative.)  This  is 

equivalent  to  saying  we  will  always  normalize  the  partial  remainder  between 

additions,  and  is  profitable  if  paths  to  perform  multiple  shifts  exist,  or  if 

the  time  for  a  given  step  is  measurably  dependent  on  whether  or  not  that  step 

requires  use  of  the  adder.   Under  these  circumstances  the  average  partial 

remainder  magnitude  is  of  interest,  since  the  average  number  of  leading  zeros  in 

R.  is  equal  to  the  average  number  of  quotient  digits  formed  per  addition  in  the 

divison  process. 

To  state  the  S-R-T  division  algorithm  formally,  the  following  defini- 
te 

tions  will  be  useful: 

a)  P.  =  the  partial  remainder  after  the  ith  step,  where 
i  =  0,  1,  . . .,  m-1 

b)  R.    =    |p.  I 

l         '    i  ' 

c)  D  =  divisor 


h 

Ref.    No.    2. 

?  Ref.    No.    10 

6 

Ref.   No.    11 

7 

Inclusion  of  end  points    is   dependent   on  choice   of  number  representation. 

1  8  T 

In  each  case,  the  fraction  portion  of  a  floating-point  number  is  implied, 
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d)   q  =  the  quotient  digit  generated  by  the  ith  step, 

with  possible  values  of  +1,  0,  or  -1 

m-1  _• 
e  )   Q  =  quot  ient  Z  2   q . 
i=0    x 

One  way  of  stating  the  recursion  relation  is  as  follows: 


l)   If  R.  <  ==,  then  P.  ,  =  2P. 
l   2'       l+l     i 


q±  =  0 


2 )   If  R .  >  ^  and 
1  —  2 


a)  the  signs  of  P.  and  D  agree,  then  P.  .,  =  2(P.  -  D) 

b)  the  signs  of  P.  and  D  disagree,  the  P.  ,  =2(P.  +D)  and 

l  i+l    v  l 

q±=  -i 

Since  the  ranges  of  —  <  P  <  1  and  —  <  D  <  1  would  imply  —  <  Q  <  2, 
corrective  initial  and/or  terminal  steps  may  be  required  to  obtain  a  quotient 
for  which  Q  is  a  proper  fraction.   It  can  be  proved  by  induction  that  the  above 
process  yields  the  correct  quotient,  since: 

P.  .  -  2(P.  -  q.D)  (2.1) 

i+l    v  l   ^i  v    ' 


and  therefore  after  m  steps  are  completed 


m-1 

P_  =  2"mP  +  D  •   Z  q.  •  2_1  (2.2) 

0      m       .  _  i 
i=0 


or 


P  =  QD  +  2_mP  (2.3) 

0  m 
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1.3  Division  as  a  Markov  Process 

The  partial  remainder  in  a  division  process  becomes  a  random  variable 
if  either  or  both  of  the  dividend  and  divisor  are  random  variables.   A  useful 
and  somewhat  surprising  relation  is  that  for  any  fixed  divisor,  the  distribution 
of  R.  is  independent  of  the  distribution  of  dividends,  for  i  suff icientely  large, 
provided  variation  in  the  dividend  density  is  bounded. 

1.3=1  Existence  and  Uniqueness  of  a  Stationary  Distribution 

In  order  to  demonstrate  that  R.  has  a  unique  stationary  distribution 

Q 
in  the  limit,  Freiman  observed  that  the  division  process  may  be  characterized 

as  a  finite,  regular  (or  aperiodic  and  irreducible)  Markov  chain.   The  states 

of  this  Markov  model  are  intervals  in  the  range  of  R. ,  each  of  which  satisfies 

i' 

the  following  requirements: 

a)  The  probability  density  function  describing  R.  is  uniform 
within  the  interval  for  all  i.  That  is,  if  xn  and  x  are 
the  end  points  of  such  an  interval, 

Pr[x  <  R  <  x  +  dx]  -  K  •  dx,    for  x  <  x  <  x  (3.l) 

where  K  is  a  constant  for  any  particular  i,  but  may  vary 
with  io 

b)  One  such  interval  in  the  distribution  of  R.  will  map 

into  an  integral  number  of  intervals  in  the  distribution 

of  R.  .  .  for  all  i.   This  means  no  new  sub intervals  are 
l+l' 

required  in  order  to  insure  that  (a)  is  satisfied  at 
each  step. 


g 

Reference  k 

10 


Proof  that  such  states  can  be  generated  is  deferred  to  Section  1.3-2 
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Denote  these  states  (S.),  j  =  1,  2,  .  ..,  N,  and  let  X(i), 

J 

i  =  0,  1,  2,  .  ..,  represent  the  stochastic  process  which  ranges  over  this  state 

space,  where  i  corresponds  to  the  step  number,  and  X(p)  =  q  denotes  the  event 

"R  is  in  state  S  ."  The  Markovian  character  of  X(i)  is  a  result  of  the  two 
P  q  . 

properties  of  all  states  listed  above.   If  X(i)  =  m  is  know  to  be  true,  the 

probability  that  X(i  +  l)  =  n,  for  example,  is  not  altered  by  information 

regarding  which  state  X(i  -  l)  assumed,  because,  due  to  the  uniformity  of  the 

density  function  within  an  interval,  this  information  would  tell  us  nothing  more 

about  the  probability  that  R.  is  in  the  subset  of  S  which  maps  into  S  .   By 

l m  n 

induction,  the  random  variable  X(i  +  l)  is  also  independent  of  all  X(p), 
p  <  (i  -  l),  as  well,  or: 

Pr[X(i  +  l)  =  n|x(i)  =  m,  X(i  -  l)  =  k,  ...,  X(o)  -  a] 

(3.2) 

=  Pr[X(i  +  l)  =  n|x(i)  =  m] 

The  rather  trivial  observation  that  the  division  algorithm  is  invariant 
from  one  step  to  the  next  enables  further  classifying  the  stochastic  process  as 
a  stationary  Markov  chain,  i.e.,  conditional  probabilities  Pr[x(i  +  l)  =  n|X(i)  =  m] 
are  transition  probabilities  which  are  independent  of  i.   The  procedure  for 
finding  the  distribution  of  R.   ,  given  that  of  R.  (or  equivalently,  finding  the 
distribution  of  the  discrete  random  variable  X(i  +  l),  given  the  distribution 
of  X(i))  can  be  formally  expressed  by  use  of  a  transition  matrix,  P,  with  dimen- 
sions NXN  for  N  states,  and  which  has  as  its  row  m,  column  n  entry  the  transition 
probability  Pr[X(i  +  l)  =  n|x(i)  =  m],  i.e.,  rows  are  identified  with  sources, 
columns  with  destinations.   Then  if  P  is  premultiplied  by  an  N-entry  row  vector, 
Jt(i),  which  describes  the  distribution  of  X(i),  i.e.,  element  Jt.  (i) 
=  Pr[X(i)  =  m],  the  product  is  Jt(i  +  l).   Since  P  is  independent  of  i,  this 
result  can  be  generalized  to: 


ir(i  +  k)  =  *(i)  •  Pk;   k  =  1,  2,  ...  (3-3) 


which  expresses  a  k-step  transition. 

The  desired  result  is  that  for  i  sufficiently  large,  fl(i)  approaches 
a  unique  final  distribution  which  is  independent  of  u(o).   If  X   denotes  this 
limit  distribution,  then 


X   =  lim(*(o)  ■  Pk)  (3.^) 

k-»o 

would  be  true  for  any  1  x  N  probability  vector,  rt(0).   Equation  (3-^)  would 
obviously  be  true  if  each  row  of  P  approached  X   as  a  limit,  since  the  sum  of 
the  elements  in  rt(0)  equals  1.   We  can  prove  that  is  precisely  what  the  limiting 
behavior  of  P  is,  and  can  state  as  a 

Theorem:   If  the  entire  interval  0  <  R.  <  1  can  be  covered  with  a  finite  number 
1 

of  intervals,  each  of  which  has  the  two  requisite  properties  of  a  state  listed 
at  the  start  of  I.3.I,  then  there  exists  a  unique  vector  X   such  that  equation 
(3»*0  is  true,  for  any  tt(0)}    each  row  of  P  approaches  X   as  a  limit. 
Proof: 

1.  The  conclusion  holds  if  the  chain  is  irreducible  and  aperiodic  (e.g.,  see 
Feller,  p.  356,  or  Kemeny,  et .  al.,  Chapt .  6). 

2.  Both  properties  are  defined  to  apply  if  there  exists  a  finite  integer  B, 

such  that  for  all  k  >  B  and  all  state  pairs  S  ,  S  ,  there  is  a  nonzero 

n7  nr 

probability  of  the  system  being  in  state  S  after  k  steps,  having  started 
in  state  S  .   (This  applies  to  the  case  n  =  m  as  well.)  Equivalently,  the 
properties  apply  if  all  elements  of  P  are  positive  for  k  >  B. 


3.   To  show  that  B  exists,  we  refer  to  equation  (2.l),  which,  if  attention  is 
confined  to  absolute  values,  is  represented  in  Figure  1.    The  proof  is 
completed  with  four  observations : 

a)  In  the  range  of  R.,  any  interval  of  length  <  —  which 
does  not  include  the  points  R.  =  .5  or  R.  =  D;  projects 
onto  an  interval  twice  its  length  in  the  range  of  R.   . 
The  same  comment  applies  to  descendents,  until  one  or 
both  of  .5  and  D  are  covered. 

b)  This  doubling  of  the  image  with  each  step  applies  also 
to  an  interval  which  includes  R.  =  .5,  provided 

R.  =  (D  -  .5)  is  not  also  included  (see  Figure  l). 
When  both  of  these  points  are  included,  clearly  D  is 
included  in  the  projected  image  on  the  range  of  R.   , 
to  which  observation  (c)  applies . 

c)  An  interval  in  the  range  of  R.  containing  D  has  a 
projected  image  expressible  as  0  <  R.   <  5,  with 
5  >  0;  on  successive  iterations  the  image  will  be 
(0,  2   •  5),  n  -   1,  2,  . ..  until  the  entire  range  of 
the  partial  remainder  magnitude,  (0,l),  is  covered. 

d)  The  first  three  observations  established  that,  having 
started  in  any  state  (nonzero  interval),  the  system  can 
reach  any  other  state,  i.e.,  irreducibility.   The  fact 
that  the  last,  single  image  of  the  sequence  listed  in 
(c)  covers  all  of  (0, l)  proves  the  existence  of  the 
integer  B,  since  this  is  equivalent  to  all  elements 
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Figure  1  derives  from  equation  (2.l),  together  with  R.  -  |P. | 


i+l 


1/2 

Ri 


FIGURE  1'    Rj+1vs.  Rj 
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of  P  being  positive  for  some  k,   Clearly  all  higher 
powers  of  P  must  also  contain  only  positive  elements., 
and  hence  the  chain  is  aperiodic. 
The  last  observation  leads  to  an  immediate 

Corollary:   All  elements  of  A,  are  positive. 

The  fact  that  the  stationary  density  function  is  nowhere  zero  in  the 
unit  interval  will  be  used  in  Chapter  2. 

1.3.2  Formulation  of  the  Chain 

First  it  will  be  shown  that  it  is  possible  to  generate  a  set  of 
intervals  covering  (0,  l)  which  meet  the  two  requirements  of  states  listed  at 
the  start  of  1.3.1;  for  a  large  class  of  dividend  distributions.  With  this 
established,  the  conditional  hypothesis  of  the  theorem  in  1„3°1  can  be  removed, 
and  for  computational  simplicity,  the  states  may  be  selected  on  the  assumption 
that  the  distribution  of  dividends  is  uniform,  since  the  final  result  is 
independent  of  the  actual  distribution  of  dividends. 

The  states  of  the  Markov  chain  for  a  given  rational  divisor,  D,  can 
be  selected  as  follows: 

a)  While  otherwise  arbitrary,  the  density  function  describing 
dividend  magnitudes  is  assumed  to  be  sufficiently  well 
behaved  that  it  can  be  approximated  in  measure  to  any 
specified  accuracy  by  a  piecewise  uniform  density,  with  all 
discontinuities  in  the  approximation  occurring  at  rational 
numbers.   (in  particular,  the  density  is  everywhere  finite; 
there  is  no  single  point  in  the  range  Rn  to  which  a  non- 
zero probability  is  assigned.  )  Denoting  the  point  set 
consisting  of  the  divisor  D  and  the  arguments  at  which 
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the  discontinuities  occur  as  -s  —  r     (each  point 
represented  in  reduced  fractional  form),  we  find  the 
least  common  denominator.   Denote  this  number  as  L. 
b)  The  intervals,  each  L   in  length,  defined  by  the 

integral  multiples  of  L  ,  are  then  the  states .   Due 
to  the  definition  of  L  in  (a),  the  density  function  for 
dividends  is  uniform  within  any  such  interval.   Discon- 
tinuity arguments  for  the  density  functions  of  R., 
i  =  1,  2t    . ..,  all  are  either  images  of  these  initial 
discontinuities,  or  of  the  points  R  =  .5  or  R  =  1. 
(See  Figure  1;  the  latter  two  are  the  points  in  the 

range  of  R.  at  which  the  R.  n  versus  R.  relation  is 
to      l  l+l         l 

discontinuous.)   Since  each  such  image  is  expressible 
in  terms  of  integral  multiples,  and  sums  and  differences, 
of  elements  in  the  set  -{- — r  s  it  is  clear  that  density 
function  uniformity  within  the  specified  intervals  is 
preserved,  and  that  the  set  of  interval  end  points 
maps  only  onto  itself. 
Other  than  showing  that  the  hypothesis  in  the  theorem  in  1.3-1  is 
satisfied,  the  above  formulation  reveals  another  interesting  property  of  the 
stationary  density  function,  viz.,  all  discontinuities  are  images  of  the  points 
.5  or  1,  since  the  uniqueness  property  implies  none  of  the  discontinuities  in 
the  R  density  and  hence  none  of  their  images  can  influence  the  limiting 
density. 
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1.3-3  Examples 

At  this  point  two  examples  may  be  found  helpful.   We  know  that 
successive  distributions  of  X(i)  are  related  by  the  equation: 

*(i  +  1)  =  K(i)  •  P  (3.5) 

for  all  i,  and  that  for  i  sufficiently  large: 

«(i)  =  *(i  +  1)  =  X  (3.6) 

Hence ^  the  solution  for  a  given  D  requires  finding  P,  then  solving  the  equation 
set 

\  =  XP  (3-7) 

To  determine  state  boundaries  we  generate  the  images  of  the  points  x  =  .5  and 

x  =  1  by  application  of  equation  (2.l)  or  Figure  1,  stopping  when  repetition 

begins.   To  determine  the  n,m  element  of  the  transition  matrix,  P,  we  answer 

the  question:   Given  an  ensemble  of  partial  remainder  uniformly  distributed  in 

state  S  ,  what  fraction  of  them  will  be  in  S  after  one  step? 
n  m 


7 

Example  1:   D  =  ^ 

a)   Images  of  the  point  R.  -  —  are: 
(i)  2(D-  |)  -| 

(ii)  2[D  -  (2D  -  l)]  =  2  -  2D  =  - 

o 
(iii)  2[2  -  2D]  =  k   -  ^D  =  2 

(iv)  2[k  -  k-D   -  D]  =  8  -  10D  =  | 

Since  the  next  application  of  equation  (2.l)  yields 

k 

q,  which  is  also  the  first  image  point  of  1,  all 

7 
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potential  locations  of  discontinuities  have  been 

determined. 

The  states  are  therefore : 


2 

? 


D 


b)   The  transition  matrix  then  is 


P  = 


.5 

.5 

0 

0 

0 

0 

0 

•  25 

•  75 

0 

0 

0 

•  5 

0 

.66... 

.33... 

0 

0 

0 

0 

l 

0 

0 

0 

c)      Solving  the  equation  set  X  =  XP,    we   obtain: 


XlA2Xl3*lAX15: 


[£  g  1  _3  _11 
L7  7  7  lU  lV 


Example  2:      D  =  ^ 


a)      Images   of    »5   are: 

(i)      2(D   -|)   -i 
(ii)     2(i')  =  | 

(in)    M§)  =  | 

(iv)     2(|  -   D)  =  | 

It 

The  first   image   of  1   is   also  — , 

5 


The   states   are 


therefore: 
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o 


D 


b) 


P  = 


•  25 


0 

0 

1 

0 

0 

1 

0 

.25 


Ull  K12  \l  \k]   =  [3  I   3  12] 


1.3.4  Limitations  of  the  Markov  Chain  Approach 

The  Markov  model  has  been  very  useful  both  conceptually  and  as  a 
vehicle  for  proving  the  existence  and  uniqueness  of  a  stationary  probability 
distribution  for  partial  remainders.   However,  when  attempting  to  use  the 
Markov  chain  to  solve  for  this  distribution,  two  significant  obstacles  are 
encountered : 

a)  The  size  of  the  linear  equation  set  to  be  solved  is 

governed  by  the  value  of  the  denominator  of  the  divisor, 
D,  rather  than  by  the  complexity  of  the  density  function 
being  sought.   This  is  a  result  of  the  definition  of 
states  which  requires  the  set  of  end  points  to  be  closed 
under  the  recursion  relation.   Thus,  in  example  1  above, 
we  find  by  normalizing  the  probability  of  each  interval 
with  respect  to  the  length  of  the  interval,  the  prob- 
ability density  function  has  only  one  discontinuity, 
although  five  equations  were  required  (see  Figure  2). 
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PR[x<R.<  x  +  dx] 


FIGURE  2:   R.  STATIONARY  DENSITY  FOR  D  =  jr 
i  9 


b)  Much  more  fundamental  is  the  fact  that  a  solution  for  a 
given  value  of  D  may  reveal  little  about  the  behavior 
of  the  density  function  in  that  neighborhood  of  the 
divisor  range.   Thus  it  happens  that  in  example  1 
if  we  change  D  to  any  other  value  in  the  interval 
(j-f    l),    the  stationary  density  still  has  only  one 
discontinuity,  as  will  be  shown  in  Chapter  2.   However, 

if  D  is  varied  even  an  infinitesimal  distance  from 

3 
D  =  —  in  example  2,  the  number  of  discontinuities  in 

the  resulting  density  function  may  become  arbitrarily 

large. 
In  Chapter  2}    the  equation  which  defines  the  steady-state  distribution 
is  examined.   By  deriving  from  it  certain  properties  which  must  characterize 
the  R.  stationary  density  functions  for  all  values  of  D,  a  method  to  systemat- 
ically generate  solutions  for  intervals  of  D,  without  the  use  of  simultaneous 
linear  equations  is  developed. 


2.   DERIVATION  OF  PARTIAL  REMAINDER  DENSITIES 

2,1  The  Defining  Equation 

The  distribution  of  partial  remainders  at  the  i+1  step  of  an  S-R-T 
division  process  may  be  expressed  in  terms  of  the  distribution  at  the  preceding 
step  as  follows : 


x 


f.(|)  +  f .  (D  +  |)  -  f*(D  -  £),    0  <  x  <  1 


fi+1(x)  =-<  1,      >  1 


0,    x  <  0 


(1.1) 


where 


)  f.(y)  =  Pr[R.  <  yJ(  for  a  given  divisor  D; 


b)  fj(y)  =< 


t±(j),        T>\ 


ft(f),  J<\ 


(1.2) 


c)  |  <  D  <  1 


This  relation  can  be  demonstrated  graphically  as  in  Figure  3>    where  the  intervals 

of  R.  which  yield  an  R.  ,  less  than  a  specified  number,  x,  are  indicated.   In 

l       J  l+l  *  '   ' 

the  limit,  if  distributions  of  successive  partial  remainders  are  equal,  let 


f(x)  ■  lim[f  .(x)] 
i— »oo 


(1.3) 


Then  we  have  equation  (l.l)  with  subscripts  deleted 
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Ri+i 


X/2 


I  /2      (D-X/2)      (O+X/2) 


Ri 


1 


FIGURE  3=    .«x)=  f(x/2)+f(D+x) 

-f(D-X/2) 
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f(|)  +  f (D  +  §)  -  f*(D  -  §),    0  <  x  <  1 


f(x)  =  "S  1,    x  >  0 


0,    x  <  0 


(1A) 


to  which  (1.2)  applies,  again  with  subscripts  deleted. 

Because  attention  is  restricted  to  piecewise  uniform  distributions, 
i.e.,  piecewise  linear  distribution  functions,  f(x)  is  differentiable  at  all 
but  a  finite  number  of  points.   Therefore,  to  find  the  density  function,  we  can 
take  the  derivative  of  both  sides  of  equation  (l.h)    (to  simplify  notation,  let 
F(x)  =  f(x)): 


where 


r 


F(x)  = 


|[F(|)  +  F(D  +§)  +  F*(D  -  §)],   0<x  <1 


=  < 


0,   otherwise 


(1.5) 


^F(y),   y>i 


F*(y)  =  < 


0, 


y       2 


Since  F(x)dx  has  the  interpretation:   Pr[x  <  R.  <  x  +  dx],  equation  (1.5)  is 
expressing  the  probability  of  a  partial  remainder  being  in  any  differential 
interval,  say  dxn,  as  the  sum  of  probabilities  that  its  predecessor  was  in  one 
of  the  (at  most  three)  mutually  exclusive  differential  intervals,  each  half  as 


long  as  dx  ,  which  map  into  dx».   Hence  we  have,  in  terms  of  a  function  whose 
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domain  is  a  continuum,   the  relation  that  was  expressed  in  terms  of  discrete 
state  probabilities  by  equation  1:(3.T)»   The  theory  of  Markov  chains  can  be 
applied  to  demonstrate  that  a  piecewise  constant  density  function  which 
satisfies  (1.5)  for  a  given  D  is  unique.   As  in  Chapter  1,  generality  will  be 
slightly  compromised  by  assuming  D  is  rational.   The  first  theorem  below  is 
instrumental  in  the  proof  of  the  second. 

Theorem  2 . la ) :   For  D  rational,  and  F(x)  a  piecewise  constant  density  function 
satisfying  (l.5)>  all  discontinuities  in  F(x)  occur  at  rational  values  of  x. 
Proof : 

1=   To  construct  a  contradiction,  assume  there  exists  an  irrational  number,  x~, 
at  which  the  conclusion  is  violated. 

2,  Since  the  number  of  discontinuities  occurring  at  irrational  arguments  must 
be  the  same  for  both  sides  of  (1.5  )j,  jumps  in  two  different  terms  on  the 
right  of  that  equation  cannot  cancel  each  other  for  x  irrational.   (As 
demonstrated  later,  such  mutual  cancellation  can  occur  for  x  rational.) 
Hence,  (1.5)  implies  that  a  second  irrational  discontinuity  point  exists 

at  one  of  2x  ,  2D  -  2xn,  or  2x„  -  2D,  depending  on  which  of  (0,  — ),  (— ,  D), 
or  (D,  l)  includes  x~. 

3.  By  induction,  the  second  such  discontinuity  implies  the  existence  of  a 
third,  etc.   The  n'th  term  in  the  sequence  is  obviously  expressible  as: 


x  =+2n  •  x^  +  C  (1.6) 

n  -      0 


where  C  is  some  rational  number  (because  D  is  rational). 
k,      Since  a  piecewise  constant  function  is  defined  to  have  at  most  a  finite 
number  of  discontinuities,  the  sequence  described  by  (1.6)  must  cycle. 


In  the  mathematical  model,  representation  of  any  partial  remainder,  rational 
or  irrational,  is  assumed  possible. 
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If  x  is  a  point  in  this  cycle,  and  recurs  after  k  steps,  we  could  infer 
from  (1.6): 


C 
x^  = 


l  +  2k 


which  is  impossible  with  x   irrational. 


Theorem  2.1h ) :   There  is  one,  and  only  one,  function  which  satisfies  equation 
(1.5)  for  a  given  rational  divisor,  D,  under  the  implied  restrictions  of  a  piece- 


wise  uniform  probability  distribution:   (i)   /  F(x)dx  =1,   (ii)  F(x)  is  piece- 
wise  constant,  and  (iii)  F(x)  >  0  for  all  x. 


Proof : 

1.  At  least  one  solution  exists,  since  the  vector  representing  state 
probabilities  of  a  Markov  chain  (i.e.,  the  solution  to  equation  1:(3.2)), 
can  alternatively  be  expressed  as  a  density  which  must  satisfy  (1.5)  at 
each  point. 

2.  Conversely,  a  solution  to  (1.5 )  for  a  given  divisor  must  be  unique  if  the 

function  can  alternatively  be  expressed  in  terms  of  state  probabilities  of 

2 
a  finite  Markov  chain,  i.e.,  if  there  exists  a  partitioning  of  the 

remainder  range  such  that  the  distribution  is  uniform  within  each  interval, 

and  the  finite  partition-point  set  is  closed  under  the  S-R-T  recursion 

relation.   Then,  for  two  alleged  solutions  to  (1.5 )  with  D  fixed,  we  could 

select  a  single  partitioning  (viz.,  the  union  of  two  partition  point  sets 

which  apply  to  the  solutions  individually)  for  both  functions,  and  thereby 

construct  a  contradiction  based  on  the  uniqueness  of  regular  Markov  chain 

state  probabilities. 


2 


While  this  uniqueness  extends  to  regular  Markov  chains  with  a  denumerably 
infinite  number  of  states,  the  a  priori  knowledge  that  the  number  of  discon- 
tinuities is  finite,  gained  by  restricting  D  to  the  set  of  rational  numbers, 
is  useful  in  the  subsequent  development  (see  Ref.  1,  Chapt.  1  for  properties 
of  chains  with  an  infinite  number  of  states). 
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3.   To  select  the  partition  point  set  hypothesized  in  2,  we  observe  that  the 
number  of  discontinuities  is  finite,  and  all  occur  at  rational  numbers 
due  to  theorem  2.1a.   Hence,  the  set  \t\}    n  =  0,  1,  . ..,    L,  where  L  is  the 
least  common  denominator  of  the  divisor  and  discontinuity  arguments,  is 
always  sufficient  for  the  partition  point  set  described  in  step  2  above, 
Q.E.D. 

The  problem  of  finding  the  stationary  density  of  partial  remainders 
is  now  reduced  to  solving  (1.5)  for  F(x),  as  a  function  of  D.   Since  no  single 
number  in  the  range  of  partial  remainders  has  a  nonzero  probability,  the  defini- 
tion of  F(x)  for  x  =  0,  1,  or  any  of  the  finite  number  of  discontinuity  arguments 
will  not  be  a  primary  concern.   Moreover,  D  will  be  restricted  to  the  open 
interval,  —  <  D  <  1,  both  for  expediency  in  derivations,  and  because  F(x)  for 
D  =  —  is  known;  to  wit:   F(x)=l,  0  <  x  <  1. 

The  definition  of  a  single  step  in  the  S-R-T  division  process  used 

3 
here  does  not  agree  with  that  used  by  Freiman.    Instead  of  defining  a  step  to 

be  either  one  binary  shift  (i.e.,  if  q.  =  0)  or  one  addition  and  shift,  Freiman 

defines  a  step  to  be  one  addition  and  partial  remainder  normalization.   If  g(x) 

denotes  the  steady-state  distribution  of  the  normalized  partial  remainder 

magnitudes,  it  is  apparent  from  what  computer  normalization  entails  that: 


'   00 

E  [f(— )  -  f(— )];    i<  x  <  1 
n       n      2 
n=0   2       2 


g(x)  =  < 


(1.7) 


x  <  I 


from  which  we  obtain: 


3  Ref .   4. 
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*(x)  =  < 


io?F(^  ^<x<1 


0,    otherwise 


(1.8) 


Hence,  information  obtained  using  one  definition  can  be  readily  compared  with 
results  for  which  the  other  was  used0   However,  two  advantages  were  found  in 
using  the  definition  adopted  here.   One  is  the  very  useful  symmetry  which 
characterizes  all  solutions  to  (l05)^  as  proved  in  Section  2=3»   The  other  is 
a  simplification  in  evaluating  the  expectation  value  for  the  number  of  shifts 
occurring  between  additions „   The  latter  point  will  be  discussed  in  Chapter  3° 

2o2  End  Point  Relations 

Two  properties  of  all  solutions  can  be  deduced  by  direct  substitution 
into  (l„5): 


F(0+)  =  2F(D) 


Ft!  -)  =  2F(1~) 


(2.1) 
(2.2) 


where;,  for  example, 


F(0+)  =   lim[F(0  +  e)] 


In  (20l),  F(D)  has  the  interpretation 


|  [F(D+)  +  F(D-)1 

if  D  is  itself  the  argument  of  a  discontinuity .   Both  (2.l)  and  (2=2)  depend 
on  the  fact  that  the  number  of  discontinuities  is  finite,  since  the  limits 


mi 
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ight  not  exist  if  0,  -,  D  or  1  were  accumlation  points  of  discontinuities 


(2.2)  also  relies  on  the  fact  that  D  =  —  is  not  being  considered, 


2„3  Symmetry 

Theorem s   All  solutions  to  equation  (1.5 )  a^e  symmetric  with  respect  to  the  point 
defined  by  coordinates  D,  F(D),  the  symmetry  extending  over  the  interval 
0  <  x  <  2D;  i.e., 

F(D  _  Q)   _  F(D)  =  F(D)  -  F(D  +  9),         0  <  9  <   D 

/l    > 
for  all  D  in  the  (— ,  1)   interval. 

Proof ; 

To  summarize  what  follows,  F  is  extended  to  a  function  with  period 

2D,  and  it  is  shown  that  the  even  portion  (excluding  an  additive  constant)  of 

the  convergent  Fourier  series  representation  of  this  extension  must  sum  to  zero 

for  all  x,  which,  due  to  Fourier  series  uniqueness,  implies  each  cosine  term 

is  identically  zero.   The  conclusion  is  an  immediate  sequel. 

1.   Define 


rF(x),    0  <  x  <  2D 


H(x)  =  <  (3.1) 


v 


F(XM0D.2D^    othe™ise 


Due  to  the  fact  that  F(x)  is  a  piecewise  constant  function  for  all  D  and  finite 
for  all  x,  F(x)  has  bounded  variation.  Therefore,  the  Dirichlet -Jordan  theorem 
applies  when  formulating  a  Fourier  series  representation  of  H(x).   If  we  define: 


h 

e.g.,  see  Zygmund;  Theorem  8.1,  Chapt .  II < 
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r2D 
bn  =  ij        sin(^)    -    F(x)dx  (3.2) 


n=±J       cos(^)   •   F(x)dx  (3.3) 


such  that  the   series   corresponding  to  H(x)   is: 

S[H(x)]    =  -|  +     S    KsinO   +  ancosOf  (3°U) 

n=l    ^  -* 

then  the  Dirichlet- Jordan  theorem  states  that  S[H(x)]  converges  uniformly 
to  H(x)  within  any  closed  interval  in  which  H(x)  has  no  discontinuities, 


and,  moreover,  if  x   is  the  argument  of  a  discontinuity 


0 


S[H(xQ)]  =  ||h(x0+)  +  H(x0-)|  (3»5) 


Therefore  the  notational  distinction  between  H(x)  and  its  series 
representation  is  unnecessary  if  we  define  H(x)  (and  F(x))  as  the  average 
of  left-  and  right-hand  limits  at  each  jump,, 
2.   Since  the  period  is  2D,  and  there  is  unit  area  under  F(x),  we  know  the 
constant  of  the  series : 


?  -  |  B  (3-6) 


To  evaluate  a  ,  n  >  0,  the  expression  for  F(x)  in  (l„5)  may  be  substituted 


into  (3.3). 
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a  = 


2D, 


cos(— ) 


1   f2D 
F(|)dx  +±J        cos(^)F(D  +  |)dx 


(3  =  7) 


2D-1 


2D, 


,nnx 


0 


cos(^)F*(D  -  f)dx 


The  limits  of  integration  in  (3° 7)  are  determined  by  the  intervals  through 
■which  x  may  vary  in  the  respective  terms  on  the  right  of  (l„5)°   The  upper 
limit  on  the  second  integral  may  he  any  value  between  2  -  2D  and  2D 
inclusive ^  all  other  limits  are  fixed „   The  flexibility  of  this  one  limit 
is  due  to  H(D  +  — )  being  zero  for  that  interval  of  x;  together  with  the 
fact  that  it  was  unnecessary  to  impose  any  nonlinearity  on  the  argument  of 
the  F(D  +  — )  term  in  (l„5).   (The  equation  defining  F(x)  could  alternatively 
be  written  as ° 


F(x)  -   |[FA(|)  +  F(D  +  §)  +  FB(D  -  |)],    |  <  D  <  1 


(3.8) 


where 


F(y),   y  <  | 


FA(y)  =  i 


0^    otherwise 


FB(y)  =  { 


F(y),   y>§ 


0y        otherwise 


no  other  restrictions  on  the  range  of  x  are  necessary, ) 

Substituting  the  series  representation  for  F  in  the  integrands  of 
(3°7);>  we  obtain  the  following  expressions: 
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.)  First  term  on  the  right  of  (3-7 ): 


p  1  °°     r 

1     /        1  /iwxv        1  v  r         /  V\   flx  / 

2D  J        2D  COs(— }   +  2  5,   lav[c°s(n  +  2}  ~  +   C03(n 

0 


2;    D  ■ 


+  bvLs-in(n  +  -)  —  +  sin(n  -  -)  -^-J  f  dx 


(3.9) 


1  .     /nit.  1  1       vi 

tt —  sin—     +y-^ar,      +  i —     Lia 
4Dnit  v  D  '        4D     2n       4it 

V=l 


sm(n  +  -)^       sin(n  -  ^  - 


V 

n  +  2 


V 

n"   2 


ri  - 


+  b 


/  Vxlt 

cos(n  +  -)- 

V 

n  +  - 


1   - 


/  VnIT 

cos(n  -  -)-• 

V 

n  -    g 


where  a  prime  on  a  summation  symbol  means  the   index,    v,    assumes   all 


positive   integer  values   excluding  2n  for  terms    in  which   (n  - 
is   a  factor.      The   same   notation   is   used  below, 
b)     Second  term  on  the  right  of   (3. 7)   is: 


Vn-1 


r  2D  oo  r 

1    /  1  /nitxv        1     v   /    -,  \VJ 

2dJ  2DCOs(—  >  +  2   ^"^r 


/  VvItX  / 

cos(n  +  ^)-^r  +  cos(n 


VvItX 

2'  D 


.      /  V  s-KX  .      ,  V  » Jtx 

bv   |  sin(n  +  -y~  +  sin(n  -  -)— 


dx 


(3-10) 


1  a  1     v.  b    fi  -  (-i)v  .  i  -  (-D 

2  a2n  "  ^  t|%    Dv  1  V  V 

V=l  ^     n  +  -  n  -  - 


n  is   constant,    V   is  the   index  of  the   series. 

6  1 

Trigonometric  identities  used  are:   cosA  cosB  =  —[cos (A  +  B)  +  cos (A  -  B)] 

cosA  sinB  =  -[sin(A  +  B)  +  sin(A  -  B)] 
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c)     Third  term  on  the   right   of    (3. 7) 


r2D-l  .  00 

1      /  ,nrtxN        1  v   /    .  O 

2D  J  cos(— )+-  Z(-l) 

0 


cos(n  +  -)-g-  +   cos(n   -  -)— 


-  b 


sin(n  +  -)-g-  +  sin(n   -   gj-jj- 


dx 


(3.11) 


sin(^)        agn(2D  -   l)         1      . 
kDnx       +  4D  +  Iwt     ~ 


fsin(n  +  g-J^       sin(n  -  ^"l 

. +  . 


n  + 


L       J 


+  b 


cos(n  +  -)- 


V 

n  +  2 


(-1)  cos(n  -  -)- 

+  ■ 


\V 


-    (-l)l 


V 

n  -   g 


Summing  the  expressions   developed   in    (3  =  9)^    (3<10)<,    and    (3«ll)j    we 
find  that    (3-7)   reduces   to: 


a     =  a^ 
n  2n 


(3.12) 


3°   As  mentioned  earlier,  the  variation  per  period  of  H(x)  is  known  to  be 
bounded.   For  variation  V  per  period  of  length  2D,,  it  can  be  shown  that 


d.     v. 

n '  —  nD 


(3.13) 


Now,  -by  induction,   (3.12 )  implies  that  a  can  be  equated  to  any  term  of 
the  form:   a,  ,  k  =  2-  °  n,  where  j  is  an  arbitrarily  large  integer .   Hence , 


(3«12)  and  (3.13)  together  imply  that  every  |a  j,  n  -  1,  2, 


is  less 


e.g.,  see  Zygmund,  Theorem  ^„12,  Chapt.  2„   He  derives  the  bound  on  a  for  a 
function  with  period  2ft,  but  the  modifications  to  his  derivation  are 
straightforward . 
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than  any  specified  positive  number.   The  question  of  summing  an  infinite 
number  of  arbitrarily  small  numbers  is  resolved  by  applying  Cauchy's  con- 
dition for  convergence,  i.e.,  given  e  >  0,  there  exists  an  integer  N  such 
that  all  partial  sums  of  the  series  after  the  Nth  are  within  e  of  the  limit 

00  00 

sum.   Because  E  a  cosf-7— )  <  E  a  ,  and   a  I  can  be  shown  to  be  less 
n    x  D   —    '  n"       '  n1 

than  e/N  by  use  of  (3.12)  and  (3.13),  it  follows  that: 


|Eancos(^)|  <  E  (3.14) 


Since  e  can  be  arbitrarily  small,  the  cosine  portion  of  the  series  must  be 
zero  for  all  x. 

For  any  integer  n,  sin(— — )  is  symmetric  with  respect  to  D,  and  there- 
for any  sum  of"  terms  of  this  form  is  also  symmetric.   Q.E.D. 


The  significance  of  being  able  to  select  2D  as  the  upper  limit  on  the 
second  integral  in  (3«7)  is  now  apparent.   If  any  other  limit  in  the  allowed 
range  of  [2  -  2D,  2D]  were  used,  the  expression  for  the  unknown,  a  ,  would  be 
in  terms  of  an  infinite  number  of  unknowns,  the  a  and  b  ,  which,  though  correct, 
would  have  presented  a  more  formidable  problem  than  equation  (3. 12). 


2.4   Corollaries  to  Symmetry 

The  symmetry  of  F(x)  with  respect  to  D  will  be  used  frequently  in  the 
remainder  of  this  chapter,  as  well  as  the  next.   The  properties  listed  in  this 
section  are  some  of  the  more  immediate  consequences  of  symmetry. 

Theorem  2.4a: 

(i)     F(D)  -i  (4.1) 

(ii)     F(x)  =|,        0<x<2D-l  (4.2) 
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(iii)  F[(2D  -  l)-]  -  F[(2D  -  l)+]  =  F(l-)  (4.3) 


Proof 


1.  The  first  assertion  follows  from  the  fact  that  all  terms  except  the 

constant  are  zero  on  the  right  of  (3°4)  at  x  =  Do  The  constant  is  given 
by  (3«6).  (Again,  the  average  of  left-  and  right-hand  limits  is  implied 
if  D  is  itself  a  discontinuity  point . ) 

2»   Since,  for  0  <  9  <   D,  symmetry  implies; 

F(D  -  0)  -  F(D)  =  F(D)  -  F(D  +  9)  (4.4) 

and  therefore,  due  to  (4.l): 


F(D  -   6)  +   F(D  +  0)  =  |,    0  <  0  <  D  (4.5) 


equation  (4.2)  follows  if  we  note  that  F(D  +  0  )  =  0  if  0  is  in  the 
interval:   1  -  D  <  9  <   D. 
3o   Equation  (4.3)  is  a  particular,  but  useful,  example  of  symmetry. 

From  properties  of  solutions  established  up  to  this  point,  we  can 

o 

deduce  the  following  about  F(x)  for  all  D  in  the  interval  r   <  D  <  1. 

a)  F(|)  -  |,  due  to  (4.2); 

b)  F(l-)  -  i=r,  due  to  (2.2)  and  (4.6a),- 

dD  (4.6) 

c)  F(D)  -  —,  due  to  (4.1) 

d)  F[(2D-  l)+]  =  ~,  due  to  (4.3)  and  (4.6t). 


8 


The  upper  limit  on  9   is  merely  to  confine  attention  to  the  interval  to  which 
symmetry  applies. 
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Because  F(x)  is  equal  to  —  at  each  end  and  the  midpoint  of  the 
interval  2D  -  1  <  x  <  1  (for  r-  <  D  <  l),  the  question  naturally  arises  as  to 
whether  or  not  solutions  for  all  D  in  (r-,  l)  are  constant  over  this  intervals 
This  question  is  one  of  many  answered  by  the 

Theorem  2„Mj;   F(x)  is  monotonic  for  x  >  0,  i.e.^  there  is  no  x  >  0  such  that 

F(xQ+)  -  F(xQ-}  >  0  (if  .7) 


Proof ' 

(A  discontinuity  described  by  (4.7)  will  be  termed  a  positive  jump> 
if  the  inequality  is  reA/ersed,,  it  is  a  negative  jump,  ) 

1.  Considered  individually^  the  respective  terms  on  the  right  of  equation  (l.5| 
would  tend  to  yield  a  positive  jump  in  F(x).,  for  x  =  x  3    if  there  were: 

a)  a  positive  jump  at  x  /2 

b)  a  positive  jump  at  D  +  Xn/23    for  x  <  2  -  2D,  or 

c)  a  negatiAre  jump  at  D  -  x r/2}    for  x_<  2D  -  1„ 

Due  to  {h,2)s    the  last  of  these  is  impossible „   Therefore,  any  positive 
jump  implies  the  existence  of  another, 

2.  Both  la)  and  lb)  above  cannot  be  true  for  the  same  value  of  xoj)  since  the 
number  of  positive  jumps  on  the  right  in  (1.5)  would  not  then  be  the  sum 
of  the  numbers  of  such  jumps  in  the  individual  terms ,    which  would  lead  to 
an  inconsistency,  since  (1,5)  is  expressing  F(x)  in  terms  of  itself. 

3.  Assume  x.  >  0  is  the  argument  of  the  largest  positive  jump  (i.e.,  the 
difference  in  (4.7)  is  maximum).   The  contradiction  is  immediate.,  since 
(1.5)  together  with  the  above  comments  require  the  existence  of  a  positive 
jump  twice  as  large.   Hence,,  there  are  none . 
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2„5  Generation  of  Solutions 

Due  to  equations  (^.2)  and  (k.6),    together  with  the  theorem  just 

proved,  we  can  confirm  that  the  stationary  partial  remainder  density  functions 

3 
for  all  divisors  in  the  range  r-  <  D  <  1  are  : 


fl 

D  ' 


0  <  x  <  2D  -  1 


F(x)  =  -J  |  D,    2D  -  1  <  x  <  1  (5.1) 

0,,     otherwise 

For  D  <  j-3   we  know  that  x  =  2D  -  1  is  again  the  argument  of  a  dis- 
continuity,, nonzero  in  size,,  since  this  jump  is  equal  to  that  at  x  -  1,  and  the 

latter  is  nonzero  since  F(x)  is  nowhere  zero  inside  the  unit  interval.   In 

x  3 

(l,3)}   this  discontinuity  appears  in  the  range  of  the  F(— )  term,,  for  D  <  r-,  and 

cannot  be  cancelled  Toy  an  equal  hut  opposite  jump  in  one  of  the  other  terms  in 
that  equation  because  F(x)  is  monotonic   Therefore,  there  will  be  a  discon- 
tinuity at  x  =  k-D  -   2}    if  2D  -  1  <  y-,  then  discontinuities  exist  at  ^+D  -  2  and 
8D  -  k;   etc. 

Similar  reasoning  applies  to  discontinuities  at  points  in  the 
(2D  -  — ,  l)  interval.   Due  to  symmetry,  a  discontinuity  in  the  range  of  either 
F(D  +  — )  or  F*(D  -  — )  is  cancelled  by  a  compensating  change  in  the  other  of 

these  terms  in  (1,5),  provided  both  terms  are  nonzero,   However,  for 

x       1 

x  >  2D  -  1,  i.e.,  T)  +  —   >  2D  -  — ,  the  F*  term  is  zero;  in  this  case  discon- 
tinuities in  the  range  of  the  F(D  +  — )  term  cannot  be  cancelled,  again  because 
F(x)  is  monotonic.   The  following  summarizes  these  observations; 

Theorem  2 , 5a :   A  discontinuity  in  F(x),  for  x  =  yn,  0  <  y     <  -}    implies  the 
existence  of  a  discontinuity  for  x  =  2yn,   Similarly,  a  discontinuity  for 
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x  =  y  ,  2D  -  —  <  y  <  1,  implies  the  existence  of  a  discontinuity  at 
x  -  2y  -  2D. 


Another  useful  conclusion  which  follows  almost  directly  from  the 
above  observations  is: 

Theorem  2.^b:   There  is  at  most  one  pair  of  discontinuities  in  F(x)  for  the 
interval  —  <  x  <  2D  -  — „ 
Proof : 

This  restriction  is  simply  a  consequence  of  the  fact  that  the  number 
of  jumps  on  both  sides  of  (1.5 )  must  be  equal  for  x  >  0.   A  pair  of  jumps  m 
the  —  <  x  <  2D  -  —  interval  cancel  each  other  on  the  right  of  (1.5),  due  to 
symmetry.   The  two  discontinuities  in  the  defining  equation  for  F(x)  (see  equa- 
tion (3°8))  are  each  the  source  of  one  negative  jump,,   and  thereby  compensate 
for  the  mutual  cancellation,  or  "loss,,"  of  two  other  jumps  on  the  right  side, 
but  no  more. 

Since  the  intersection  of  the  domains  of  the  three  terms  on  the  right 
of  equation  (3=8)  is  null,  while  the  union  covers  the  unit  interval  (excluding 
the  points  x  =  D,  x  -  —)f    each  discontinuity  in  F(x)  appears  in  precisely  one 
of  these  individual  terms.   Due  to  the  two  additional  discontinuities  intro- 
duced by  the  definitions  of  F„  and  F^  in  (3.8),  the  number  of  discontinuities 
in  F'(x),  x  >  0  is,  in  general,,  two  less  than  the  sum  of  the  numbers  of  discon- 
tinuities in  these  individual  terms.   This  difference  can  be  accounted  for 
either  by  the  mutual  cancellation  of  jumps  in  the  (— ,  2D  -  — .)  interval, 
described  in  theorem  2.5b,  or  by  two  pairs  of  jumps  in  the  individual  terms 


r) 


viz. :   jumps  at  x  -  2D  -  1  and  x  =  1  due  to  the  zero-going  of  F(D  -  — )  and 
F(-),  respectively. 
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summing  to  a  single  pair  in  F(x).   The  stationary  density  function  for  D  =  -jr 
is  an  example  of  the  latter  (see  Figure  k).      Since  discontinuity  arguments  are 
all  linear  expressions  in  D  (e.g..,  2D  -  1,  h  -   GV)}    etc.),  this  effective  coin- 
cidence of  two  pairs  of  discontinuity  arguments  in  general  occurs  only  for 
isolated  values  of  D,  which  will  be  called  the  "exceptional  divisors,"  denoted 

[D„}„   To  be  precise,  the  set  fD^l  is  defined  to  include  all  divisors  to  which 
v  E  E 

(5. 2d)  below  applies. 

We  are  now  in  a  position  to  find  the  arguments  of  all  discontinuities 
in  a  very  straightforward  way,   A  sequence  {x.}  is  formed.,  using  theorems  2.5a 
and  2.5b,  starting  with  x  -  2D  -  1: 

a)   If  x.  <  — ,  then  x.  ,  =  2x.  . 
l   2.3  l+l     l 


(5.2) 


b )  If  x .  >  2D  -  \,    then  x .  .  -  2D  -  2x . . 

l        2J  l+l  l 

c)  If  —  <  x.  <  2D  -  — ,  the  sequence  is  complete, 

2—i        2J 

d)  The  sequence  is  otherwise  terminated  at  the  ith  step 

if  there  exists  an  x . ,  j  <  i„  such  that  x„  =  x. . 

J  J    i 

Elements  of  a  set  fx!},  where  x!  =  2D  -  x. ,  are  then  also  discon- 

±JJ  l         i^ 

tinuity  arguments,  because  of  symmetry.   The  fact  that  [x.)    and  {x!}  include  all 
discontinuity  arguments  is  evident  from  considerations  similar  to  those  used 
to  prove  no  jumps  are  positive:   in  general,  any  discontinuity  other  than  the 
pair  resulting  from  discontinuities  in  the  defining  equation  (at  x  =  1; 
x  =  2D  -  l)  implies  the  existence  of  a  predecessor  (in  the  sense  of  the  sequence 
defined  by  (5.2))  at  which  the  jump  is  twice  as  large,  because  of  the  coefficient, 


10      1 

xi  ~  p  anc*  xi  ~  ^  are  limiting  cases  where  two  otherwise  distinct 
discontinuities  become  coincident. 
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F(X) 


5/3 
5/4 
5/6 
5/12 

2 

r 

4                   .8         1 

FIGURE  4q:   SOLUTION  FOR  D  =  3/5 


1/2  F(X/2) 
(X<1) 


1/2  F(D+X/2) 


FIGURE  4b 


FIGURE  4c 


1/2  F*(D-X/2) 


FIGURE  4d 
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-,  on  the  right  in  (1.5).    Since  the  two  just  cited  are  the  only  x.  which  do 
not  require  such  a  predecessor,  all  discontinuities  must  he  descendent  from 
these,  in  the  sense  of  (5-2) 

With  the  number  of  discontinuities  known,  the  size  of  the  respective 
jumps  can  also  be  easily  evaluated  if  we  ignore  the  set  of  divisors,  (-D^}. 
Because  F(x)  is  monotonic,  the  sum  of  the  jumps  for  x  >  0  must  be  —  (see  (k.2j). 
Moreover,  the  sizes  of  the  jumps  in  the  sequence  {x.}  form  a  geometric  progres- 
sion, i.e.,  the  jump  at  x  =  x .  -.is  twice  as  large  as  the  jump  at  x  =  jc.    because 
of  the  coefficient  of  —  in  equation  (1.5 )  -   Denote 

a.  =  F(x.-)  -  F(x.+)  (5.3) 


then  for  T  discontinuity  pairs: 


£  =  2  Lf2  Z  a^)1  =  1^(1  -  (|)T)  (5J0 

i=l       i=0 


or 


T-2 
2 


1    (2T  -  1) 


(5.5) 


(The  coefficient,  2,  in  (5-^)  accounts  for  the  fact  that  the  {a.}  are  half  of 
the  jumps.)  Therefore,  the  jump  at  x .  is: 

a,  =  (I)1'1  ■  a  =  2T"X"  (5=6) 

1    2        1       D(2T  -  1) 


11 


If  the  divisor  is  in  {D  },  or  either  of  the  cases  x.  =  D  or  x-  =  p  occur,  the 
predecessor  requirement  applies  to  all  except  one  "generation  of  discon- 
tinuities, the  exception  resulting  from  the  fact  that  one  jump  is  half  the 
sum  of  two  others. 
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jumps  at  corresponding  symmetric  discontinuities  are,,  of  course, 
equal. 

The  solutions  derived  using  (5.2)  and  (5° 6)  apply  to  intervals  of  the 
divisor  range.,  rather  man  individual  values  of  the  divisor  (excluding  solutions 
for  [B  } )}    since  there  are  T  pairs  of  discontinuities  in  F(x)  for  all  divisors 
which  result  in  x  being  included  in  the  interval  {—s    2D  -  — )0   As  demonstrated 
below.,  there  are  an  infinite  number  of  such  intervals-,  eacn  of  which  involves 
a  distinct  sequence  of  operations  when  applying  (^.,2)y    though  several  non- 
adjacent  intervals  may  result  In  the  same  number  of  discontinuities „   Equation 
(5«l)  gives  the  solution  for  one  of  these  intervals,, 


To  formulate  systematically  the  R_.  distributions  for  —  <  D  <  f^ 


1  ^  T  .3 

*  2  <  D  <  Tp 


consider  the  intervals  of  D  for  which  (— )    <  (2D  -  l)  <  (— )  9    n  =  1;  29    . ..., 

x>  ts   i    -     1   /lvn+2  .  ^  -  1    ,l\n+l 
is  true;  denote  tnese  intervals  of  D  (viz. :  —  +  {—}  <  D  <  —  +   {^)        3 

n  =  1,  2f    . . . )  as  I  ,  n  =  I,  2}  .  ...j  respectively.   To  express  a  sequence  of 

discontinuities   generated  by  use  of  (5-2)^  the  following  definitions  will  be 

useful : 

N("k.)  =  k  consecutive  normalizing  operations  (i.e.,  (5.2a)) 

S  =  one  application  of  (5 -2b) 

(Listing  tne  operations  required  to  form,  the  discontinuity  arguments^ 
rather  than  listing  the  arguments  themselves s   will  enable  observing  similarities 
in  different  intervals  of  D  which  can  then  be  used  as  a  basis  to  form  infinite 
families  of  solution  intervals.) 

Then  every  1   includes  a  sublnterval  for  which  the  sequence  of  opera- 


tions in  (5.2)  is  N(n)j,  i.e.,  the  first  x.  >  -  also  satisfies  (5. 2c).   Tn 


general,  I  will  also  include  a  subinterval  for  which  the  discontinuity  locations 
n 
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in  F(x)  are  determined  in  (5.2)  by  the  sequence  N(n)S;  another  for  which  the 
sequence  is  N(n)SN;  etc.   Table  1  is  a  list  of  all  subintervals  of  I  ,  n  >  2} 
for  which  the  divisor  results  in  n  +  5  discontinuity  pairs  or  less  (excluding 

{D  ]).   The  following  two  examples  are  provided  to  show  how  Table  1  was 

E 

derived;  the  solution  for  I  appears  in  the  first  example 

Example  1,   Consider  the  sub interval  of  I  for  which  the  sequence  in  (5°2)  is 
N(n);  this  means  the  only  discontinuities  in  F(x),  for  all  values  of  D  in  this 
subinterval  are  at:  x  =  2D  -  1;  x  =  2 (2D  -  l);  ...;  x  =  2n(2D  -  l), 
together  with  the  respective  2D  complements.   For  this  portion  of  I  ,  the  fol- 
lowing is  true : 


|  <  2n(2D  -  1)  <  2D  -  \  (5.7) 


or 


n+1  n+1 

J-±  <  D  <      "    ,        n  =  1,  2,  ...  (5.8) 

n+2   -      n+2   ^  >  > 


(See  Figure  5«)  Note  that  I  ,  .625  <  D  <  .75>  is  covered  completely  in  this 
family  of  solutions.   Since  the  number  of  discontinuity  pairs,  T,  is  two  for 
I  (in  general,  T  exceeds  the  number  of  steps  in  (5°2)  by  one)?  we  find  by  use 
of  (5.6); 
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TABLE  1.   SOLUTION  INTERVALS  IN  THE  DIVISOR  RANGE 


2r(5o2) 

Interval 
Bounds 

Examples 

(to  eight  si 

gnificant  fi 

gures ) 

Sequence 

n  =  2 

n  =  3 

n  =  4 

n  =  5 

1)   N(n) 

2n+1  -  1 
2n+1  -  4 

.5833.... 

o.53571^8 

.516..... 

. 50806452 

2n+1  +  1 
n+2 

»5625 

.53125 

.515625 

.5078125 

2)   N(n)S 

2n+2  -  1 
2n+3  _  8 

0625 

o55357l43 

.525 

.51209677 

2n+2  +  1 
2n+3  -  4 

060714286 

°55 

.52419355 

. 51190476 

3a)  N(n)SN 

2n+3  -  1 
2n+k    -  12 

0 59613846 

.54310345 

.52049180 

.51 

2n+3  -  1 
2n+k  -   8 

. 58928571 

o54l6.o.o 

. 52016129 

.50992603 

3b)  N(n)SS 

2n+3  -  1 
2n+k  -   16 

.5625 

.52916.00 

0  51411290 

2n+3  +  1 
2n+k   -  12 

......... 

. 56034483 

0  52868852 

.514 

4a)  N(n)SNN 

2U+k    -    1 
2n+5  -  20 

......... 

.52813559 

.51829263 

0 50896414 

2n+k   +  1 
2n+5  -  16 

.5375 

0 51814516 

.50892857 

4b)  NfnjSNS 

2n+k   -    1 
2n+5  -  2k 

060576923 

o5474l379 

.52254098 

.511 

2n+k   +  1 
2n+5  -  20 

060185185 

0 54661017 

.52235772 

0  51095618 

4c)  N(n)SSN 

2U+k   -  1 
2n+5  -  28 

.55701754 

. 52685950 

. 51305221 

2n+k   +  1 
2n+5  .  24 

000000.00 

.55603448 

.52663934 

.513 

4d)   N(n)SSS 

2n+4  -  1 

2n+5  -  32 

0  0  0 ....  0 . 

......... 

.53125 

.51512097 

2n+4  +  1 
2n+5  -  28 

......... 

.53099174 

.51506024 
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TABLE  1.   SOLUTION  INTERVALS  IN  THE  DIVISOR  RANGE  (CONT'D) 


2:(5»2) 

Interval 
Bounds 

Examples 

(to  eight  si 

gnificant  fj 

gures  ) 

Sequence 

n  =  2 

n  =  3 

n  =  4 

n  =  5 

5a)  N(n)SNNN 

2n+5  -  1 

2n+6  -  36 

.51720648 

1  . 50844930 

2n+5  +  1 
2n+6  -  32 

......... 

.51713710 

.50843254 
i 

5b)   N(n)SNNS 

2n+5  -  1 
2n+6  -  ^0 

.  58796296 

. 54025424 

. 51930894 

.50946215 

2n+5  +  1 

2n+6  -  36 

. 58636364 

o53991597 

.51923077 

o 5094433^ 

5c)   N(n)SNSN 

2n+5  -  1 
2n+6  -  44 

059905660 

. 54487179 

.52142857 

;  o 51047904 
1 

2n+5  +  1 
2n+6  -  40 

.5972,. ,o 

„ 54449153 

.52134146 

j  .51045817 

t 

5d)  N(n)SNSS 

2n+5  -  1 
2n+6  -  U8 

. 5495689T 

o52356557 

|  ^5115 

2n+5  +  1 
2n+6  -  44 

.5491^-530 

.  52346939 

.51147705 

5e)  N(n)SSNN 

2n+5  -  1 
2n+6  -  52 

„  55434783 

„525720l6 

1   .51252505 

1 

2n+5  +  1 
2n+6  -  48 

......... 

o55387931 

.52561475 

.5125 

5f)   N(n)SSNS 

2n+5  -  1 
2n+6  -  56 

.55921053 

.52789256 

,51355422 

2n+5  +  1 
2n+6  -  52 

o55869565 

.5277.o,. 

. 51352705 

5g)   N(n)SSSN 

2n+5  -  1 
2n+6  -  60 

.53008299 

.51^58752 

2n+5  +  1 
2n+6  -  56 

1 

.52995868 

.51^55823 

5h)   N(n)SSSS 

Qn+5 
2_-_l_ 

2n+6  -  64 

— -■■  -\ 

,515625 

2n+5  +  1 
2n+6  -  60 

......... 

,51559356 

-1+0- 

lO 

• 

o 

CM 
1 

q 

cvl 

V 

N. 

Q 

V  > 

lO 

Q 
CVJ 

ID 

• 

— 

c 
CVJ 

lO 
CM 

(& 

\ 

^ 

CM 

^ 

(0 

^ogf 

\ 

CO 
CM 

to 

9T=gg^ 

""       -   ■X-cz: 

ro 
lO 

CM 

ro 

V, 

ro 
ro 


0> 


CO 


CD 


to 


-la- 


r . 


1 

D  > 


0  <  x  <  2D 


—  ,    2D  -  1  <  x  <  min(lj-D  -  2;  2  -  2D) 


F(x)  ==  ^  —  ',        min(^D  -   2}    2  -  2D)  <  x  <  max(^D  -  2,  2  -  2D)   (5.9) 


7==:  ,   max(>D  -  2,  2  -  2D)  <  x  <  1 


\_ 


0,,     otherwise 


For  D  in  any  given  solution  subinterval,,  x  and  x'  =  2D  -  x   (and  only  these  two 
discontinuity  arguments)  may  be  related  by  any  of:   x  >  x',  x'  >  x  _,  or 
x'  =  x  =  Do   This  is  the  reason  for  the  "max"  and  "min"  in  (5°9)° 


Example  2  „   If  the  sequence  in  (5-2)  is  N(n)SS.,  the  following  relation  is 
satisfied  by  D° 


~  <  2{2n(2D  -  l)  -  D]  -  D)  <  2D  -  | 


(5«io) 


or 


2n+3  +  1        2n+3  -  1 
- — , —  <  D  - 

2U+k   -  12 


2n+k   -  16 


(5=11) 


There  are  two  reasons  why  any  given  sequence  of  operations  in  (5-2) 
may  not  correspond  to  a  solution  interval  in  the  range  of  D. 

a)  Some  intermediate  term  may  necessarily  satisfy 

—  <  x.  <  2D  -  — }    if  the  implied  final  term  also  does, 

b)  The  sequence  may  not  be  self-consistent;  e.g.,,  if  a 
sequence  begins  with  N(k)S  . .,,    then  (2D  -  l)  >  (— ) 


is  true,  and  therefore  an  inconsistency  would  exist 
if  k  +  1  (or  more)  successive  normalizations  occurred 
at  any  point  subsequently,  because  2D  -  1  is  the 
smallest  discontinuity  argument „ 
In  the  example  being  considered,,  N(l)SS  cannot  correspond  to  a  solution 
because  all  values  of  D  for  which  only  one  step  is  required  to  normalize  2D  -  1 
are  part  of  the  interval  in  which  N(l)  is  the  complete  sequence „   N(2)SS  cannot 
yield  a  solution  since,  if  the  first  three  steps  are  N(2)S,  it  can  easily  be 
shown  that  x,  <  2D  -  —3    and  therefore  the  sequence  must  either  terminate  after 
N(2)S,  or  be  followed  by  a  normalization  as  the  next  step. 

For  all  n  >  3^  equation  (5.1l)  correctly  defines  intervals  of  D  for 
which  the  discontinuities  in  F(x)  are  generated  by  N(n)SS  in  (5»2),  since 
neither  of  the  two  reasons  for  exclusion  are  violated  in  these  areas , 

To  facilitate  deciding  whether  a  given  sequence  of  operations  when 
applying  (5-2)  corresponds  to  some  interval  of  D,  there  are  four  simple  rules 
to  test  for  consistency. 

1)  In  a  sequence  starting  with  N(k)S,  the  maximum  number  of  successive 
normalizations  which  can  occur  anywhere  in  the  sequence  is  k.   This  is 
because  2D  -  1  is  the  smallest  discontinuity  argument 0 

2)  A  sequence  starting  with  N(k)S  cannot  also  terminate  with  N(k).   Consider 

Figure  5  for  a  given  value  of  D.   If  2  (2D  -  l)  >  2D  -  -,  then  any  other 

x.  occurring  in  the  (-)    <  x  <  (— )   interval  will,  after  normalization, 

also  exceed  2D  -  — .  again  because  all  x.  >  2D  -  1„ 
2'  1  — 

The  other  two  restrictions  are  simply  duals  of  the  above „   Recalling 
that  (5-2)  is  generating  half  of  the  total  set  of  discontinuity  arguments,  it 
may  be  noted  the  sequence  of  operations  to  generate  the  2D  complement  set, 
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{x! ),  starting  with  x!  =  1,  is  the  dual  (i.e.,  interchange  N's  and  S's)  of  the 

sequence  to  generate  (x.}.   To  determine  the  bound  on  the  number  of  successive 

S's  which  may  occur  for  D  in  a  given  interval,  I  ,  we  assume  D  is  such  that 

some  x.  is  almost  equal  to  one.   Rather  than  examining  bounds  on  x .  .,  .  etc.. 
$  J+l      ' 

attention  may,  at  that  point,  be  shifted  to  the  2D  complement  sequence  x! , 

J 

x'.  ,  i  etc.,  and  the  limits  on  the  number  of  successive  normalizations  in  the 
J+l 

complement  sequence  determined.   From  these  considerations  one  obtains: 
3)  A  sequence  starting  with  N(k)S  can  have  at  most  k  successive  S's. 
k)     A  sequence  starting  with  N(k)S  cannot  terminate  with  a  subsequence  of  k 
successive  S's. 

After  rejecting  by  inspection  those  sequences  which  are  not  consistent, 
the  final  test  for  the  validity  of  a  proposed  sequence  of  operations  in  (5»2) 
is  comparison  of  the  interval  of  D  implied  with  intervals  corresponding  to  all 
subsequences  (formed  by  truncation)  to  determine  whether  an  intermediate  x. 
was  in  the  (— ,  2D  -  — )  interval  as  well  as  the  implied  final  one. 

It  is  therefore  only  necessary  to  find  the  smallest  value  of  n,  say 
n  ,  for  which  a  given  sequence  in  Table  1  correctly  describes  a  solution 
interval,  and  the  expression  then  applies  for  all  n  >  n  .   Solutions  in  Table  1 
are  arranged  in  the  order  of  increasing  complexity,  such  that  all  subintervals 
of  I  for  which  F(x)  has  k  discontinuity  pairs  are  listed  before  those  with 
k  +  1,  Moreover,  by  ordering  those  which  have  the  same  number  of  discontinuities 
for  I  according  to  a  binary  weighting,  such  that  S  corresponds  to  1  and  N  to 
zero,  it  happens  that  constants  in  the  denominators  of  successive  upper  bounds 
(as  well  as  successive  lower  bounds)  form  a  readily  discernible  arithmetic 
progression.   If  desired,  this  progression  may  be  used  to  extend  the  table 
indefinitely. 
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2.5.1  Properties  of  Solution  Intervals 

Though  by  no  means  unique,  the  ordering  of  solution  intervals  given 
by  Table  1  displays  one  of  several  interesting  characteristics  of  the  solutions 
Aside  from  the  finite  number  of  I  for  which  a  given  form  of  sequence  does  not 
apply,,  all  I  include : 

a)  one  subinterval  which  yields  n  +  1  discontinuity  pairs 
in  F(x),° 

b)  one  for  which  there  are  n  +  2  pairs.; 

c)  two  for  which  there  are  n  +  3  pairs ; 

d)  four  for  which  there  are  n  +  k   pairs; 

and.  in  general,, 

e )  2  for  which  there  are  n  +  j  +  2  pairs . 

I,  represented  by  the  column  for  n  =  5j  is  an  example  of  this,,  up  to  ten 
discontinuity  pairs „ 

A  second  property  which  is  apparent  from  Table  1  is  that  solution 
intervals  are  of  the  form 


<D<—T^-i— -  (5.12) 


2m+1  -  B       2m+1  -  (B  +  k) 


where  B  and  m  are  integers „   This  can  be  proved  by  observing  that  for  an 
arbitrary  sequence  of  operations  in  (5*2),,  the  inequalities  which  determine 
the  solution  interval  (e.g.,  (5-07),  (5.10))  are  of  the  form: 

|  <  2'3(2D  -  1)  -  kD  <  2D  -  |  (5-13) 

where  j  and  k  are  integers.   (5.13)  then  yields  an  expression  of  the  form  of 
(5.12). 
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Consider  now  the  problem  of  finding  contiguous  solution  intervals, 
rather  than  the  successively  complex  but  usually  nonadjacent  intervals  of 
Table  1.   Let  d  <  D  <  d  represent  a  solution  interval  for  which  the  discon- 
tinuity arguments  generated  by  (5.2)  are  y  3    y  3    „„„,,  y  3    each  expressed  as  a 
function  of  D„   For  D  in  the  region  immediately  to  the  right  of  this  interval., 
y  exceeds  2D  -  — 3    by  definition  of  an  interval  boundary  (see  (5°2c))„   The  only 
difference  in  the  first  T  discontinuity  arguments  for  D  =  d  +,    as  compared  with 
those  for  d  <  D  <  d  3   would  be  if  some  y.3    j   <   T,  satisfied  —  <  y.  <  2D  -  —  3 
resulting  in  an  interval  of  fewer  discontinuitites .   Otherwise  stated,  since 
all  y.  are  of  the  form  bD  -  C,  where  b  and  c  are  positive.,  the  only  change  in 
the  relations  of  the  y,  to  the  thresholds  in  (5.2a)  and  (5.2b)  which  can  occur 
by  increasing  D  is  that  one  (or  more)  of  the  y.  cross  — .   If,,  however,  such  a 
y„,  j  <  T  does  not  exceed  —  3    then  d  is  an  accumulation  point  of  solution 
intervals!   To  prove  this5  we  note  that  for  D  =  d  +  £3    and  £  small  but  positive: 


yT  =  (2D  -  |) +  p  •  2Te 


where  —  <  p  <  1,,  and  therefore 


T+1 
yT+1  =  2[yT  -  D]  =  2D  -  1  +  p  •  2   e 


For  €  sufficiently  small,  the  process  defined  in  (5»2)  can  be  made  to  approxi- 
mate an  infinite  cycling  for  an  arbitrary  number  of  iterations  since  y 
approximates  y  3    and  the  perturbation  introduced  by  e  causes  the  term  which 
corresponds  to  y  in  each  iteration  to  move  further  from  the  2D  -  —  threshold „ 
To  summarize,  the  above  considerations  reveal  any  solution  interval 
is  bounded  on  the  right  by  either: 

a)  an  accumulation  point  of  solution  intervals,  or 

b)  a  solution  interval  yielding  fewer  discontinuities. 
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There  is  never  a  single  adjacent  interval  on  the  right  which  yields  a  greater 
number  of  discontinuitites . 

Next  consider  the  relation  of  solution  interval  d  <  D  <  d  to  the 
adjacent  interval  on  the  left,  again  denoting  the  terms  generated  by  (5.2) 
for  d  <  D  <  d p  as  y  ,  y  ,  . ..,  y  .   Due  to  the  definition  of  an  interval 
boundary,  as  D  is  varied  from  d  +  to  d  -,  the  y  term  moves  out  of  the  (— 3 
2D  -  — )  interval;  moreover,  none  of  the  y.,  j  <  T  replace  y  in  that  interval 
since  this  would  imply  that  d  <  D  <  d   is  bordered  on  the  left  with  an  interval 
of  fewer  discontinuity  pairs,  which  contradicts  what  was  just  established  above . 
Therefore,  the  inequality  relation  of  each  y.,  j  <  T,  with  respect  to  the 
(— ,  2D  -  _-)  interval  is  not  altered  as  D  changes  from  d  +  to  d,  -,  which  means 

the  first  T  terms  generated  in  (5°2)  for  D-d  -  are  the  same  as  those  for 

1        T 
d  <  D  <  d  „   If  D  =  d   -  G,  with  G  small,  then  y  =  —  -  p  •  2   •  €,  where 

|  <  p  <  1;  yT+1  =  1  -  P2T  •  £}    yT+2  =  2  -  2D  -  p  '  2T+1£;  . . . .   Starting  with 

y   ,  these  terms  may  be  expressed  relative  to  the  2D  complements  of  the  first 

T  T+2 

T  terms  as  follows:   yT+1  =  y|  -  p  •  2  e,  yT+2  =  y^   -   p  •  2   e,  . .., 


'prn  ""  </i 


?T  T  11 

'  -  p  •  2   Go   Now,  since  y  is  p  °  2  G  away  from  the  (— ,  2D  -  — )  interval 


y'  is  as  well,  due  to  symmetry,  and  therefore  these  expressions  reveal  that 
y   must  be  inside  the  interval.   Therefore,  d  <  D  <  d  ,  the  generic  solution 
interval,  is  bounded  on  the  left  by  an  interval  with  twice  the  number  of  dis- 
continuity pairs . 

By  induction,  the  results  of  the  last  three  paragraphs  yield  the 

Theorem  2.5c:   Every  solution  interval  is  an  element  of  an  infinite  sequence 
adjacent  intervals,  each  of  which  results  in  half  the  numbers  of  discontinuity 

pairs  in  F(x)  as  the  interval  at  its  left.   Except  for  the  sequence  beginning 

3 
with  the  f  <   D  <  1  interval,  the  first  interval  in  each  such  sequence  is 

bordered  on  the  right  by  an  accumulation  point  of  solution  intervals. 
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3 
An  example  of  such  an  infinite  sequence  begins  with  ■?  <   D  <  1;  for 

which  F(x)  has  one  discontinuity  pair.   The  next  four  intervals  are: 

a)  §<  D<  j2 

and 

d)  108124  <  D  <  ilil20^in  WhlCh  there  are  2j    k>    &J    and  l6 
discontinuity  pairs,  respectively. 

The  discontinuity  arguments  which  satisfy  the  —  <  x .  <  2D  -  -  relation  for  the 
last  three  examples  are  shown  in  Figure  6,    which  illustrates  the  rapid  con- 
vergence of  such  a  sequence,,   It  happens  that  all  the  solution  intervals 
appearing  in  Table  1,  except  those  described  by  N(2)Si  N(3)SS^  N(^)SSS.  and 
N(5)SSSS,  are  examples  of  first  intervals  in  Infinite  sequences.  Due  to  tne 
above  theorem,  it  is  only  necessary  to  observe  whether  any  of  the  intervals 
resulting  in  fewer  discontinuitites  than  the  interval  in  question  is  adjacent, 
to  decide  whether  a  given  interval  is  first  in  such  a  sequence . 

Parenthetically,  it  may  be  noted  that  equation  (5ol2)  provides  a  very 
straightforward  method  for  computing  boundaries  on  successive  adjacent  intervals 
in  these  infinite  sequences .   As  an  example,  -^k   =  — ^o —  is  of  the  form  of  a 

lower  bound  in  (5.12).   To  express  the  same  number  in  the  form  of  an  upper 

17   2^  -  1   2^-1 

bound.  7-5  °   — j —  =  — rT^r~°      Then  the  lower  bound  on  the  same  interval  is 

20    2^"  -  1      42U 

25^  +  2 

Tpk — ~T}    again  using  (5. 12).   In  addition,  the  list  of  operations  when  applying 

(5«2)  for  each  interval  of  D  in  such  a  sequence  may  be  generated  by  application 

of  the  following  rule : 
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If  u(N,S)  denotes  the  N-S  sequence  corresponding  to  any 

given  solution  interval,  and  u(N, S)  is  defined  to  be  the 

same  sequence  with  N's  and  S's  interchanged,  then  the  (5-l4) 

solution  interval  to  the  left  of  the  one  in  question  is 

characterized  by  u(N, S)Nu(N, S). 

This  rule  is  merely  a  restatement  of  an  observation  which  led  to  Theorem  2.5c, 

v 

where  it  was  noted  that  with  K  discontinuity  pairs,  the  second  —  terms 
generated  by  use  of  (5.2)  closely  approximate  the  2D  complements  of  the  first 

V 

— .  Applying  it  to  the  four  intervals  mentioned  above,  one  obtains: 

|<D<f:  N 


|J  <D  <|:  NNS 


fir  <  D  <  W      NNS  N  SSN 


65337  <Ti<    65535 
108124      108120 


As  a  final  observation,  an  immediate  corollary  to  Theorem  2.5c  is 
that  all  solution  intervals  with  odd  numbers  of  discontinuity  pairs  have  an 
accumulation  point  of  discontinuitites  as  an  upper  bound,  i.e.,  such  intervals 
are  first  in  an  infinite  sequence  of  the  type  described  by  Theorem  2.5c 


2.5.2  Comments  on  the  Set  (D_,} 
t  — 

The  exceptional  divisors,  {£>„},  defined  to  be  those  values  of  D  for 
which  (5.2d)  is  satisfied,  may  be  recognized  as  one  of  three  distinctly  different 
types  of  limit  points  (or  accumulation  points)  which  would  exist  if  Table  1  were 
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extended  indefinitely.   One  type  of  limit  point  is  the  right  edge  of  certain 
solution  intervals;  for  D  arbitrarily  close  but  greater  than  such  a  boundary, 
the  sequence  of  operations  in  (5.2)  appears  to  cycle  as  discussed  in  the 
previous  section.   However,  the  {D  }  are  accumulation  points  of  solution 
intervals  on  either  side,  unlike  these  interval  boundary  points.   The  third  type 
of  limit  point  is  that  to  which  an  infinite  sequence  of  contiguous  intervals 
converges  (see  Theorem  2.5c).   This  type  of  accumulation  point  of  solution 
intervals  is  distinguished  from  the  other  two  by  the  fact  that  the  sequence  of 
operations  when  applying  (5  =  2)  does  not  cycle.   See  (5.1^-)  and  the  examples  which 
follow  it. 

To  generate  a  few  subsets  of  [D },  consider  the  following  examples 
of  cycling  which  could  arise  when  applying  (5.2): 

a)  x  =  2D  -  1;  x  -  ^D  -  2;  x  -  8D  -  k; 
x,  =  ikl)   -  8  =  kD  -  2,    i.e.,  D  =  |. 

b)  x  =  2D  -  1;  x  =  kD  -   2;  x  =  8D  -  k}   x^  =  l6D  -  8; 

x  -  30D  -   l6  =  kO  -  2,    i.e.,  D  =  yL 

The  first  of  these  is  an  example  of  the  subset  of  {D  }  for  which  the  sequence 

ill 

of  operations  is  N(n)SNSN  ...;  n  =  2,  3,  k,  ....  The  second  is  an  example  of 
N(n)SMSM  . ,.}  n  =  3,  k,  ....  The  values  of  D  in  the  first  of  these  subsets 
satisfy  the  equation 

2n(2D  -  1)  =  2[2n+1(2D  -  l)  -  D],    n  =  1,  2,  .  .  .        (5.1^) 
and  therefore  the  subset  is : 


D-  3-^— ,       n-1,2,    ...  (5.0.}) 
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The  values  of  D  in  the  second  subset  are  described  by: 

2n(2D  -  1)  =  2[2n+2(2D  -  l)  .-  D]  (5.1$) 


or 


D  = Li^J ,    n  =  1,  2,  ...  (5.17) 

7  •  2    -  2 


Both  (5.16)  and  (51-7)  can  be  shown  to  be  proper  subsets  of  the  set  in  {D  } 


2n  -  1 


where  n  assumes  all  integer  values  excluding  powers  of  two. 
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3.   SHIFT  AVERAGE 

3.1  Equation  for  <S> 

While  the  density  of  partial  remainders,  as  a  function  of  D,  is 
interesting  in  itself,  its  practical  value  lies  in  evaluating  what  fraction  of 
the  total  number  of  quotient  digits  produced  are  nonzero.   This  figure  of 
merit  is  frequently  expressed  in  reciprocal  form  as  the  shift  average,  i.e., 
the  average  number  of  binary  shifts  in  the  division  process  between  successive 
uses  (addition  or  subtraction)  of  the  divisor.   The  shift  occurring  as  part  of 
a  step  producing  a  nonzero  quotient  digit  is  counted  when  computing  the  shift 
average .   Define : 

S(D)  =  S  =  random  variable  denoting  the  shift  count  between 
successive  uses  of  divisor  D; 

<S>  =  expectation  value  of  S. 

Since  the  distribution  of  R  is  independent  of  the  dividend  (for  i 
sufficiently  large),  the  distribution  of  S  is  also.   This  independence  of  R. 
was  based  on  the  assumption  that  the  probability  density  describing  dividends 
is  piecewise  constant,  and  finite  (Section  1.3.2).   Therefore,  for  a  divisor 
D  ,  <S(Dn)>  =  C  means  that  the  shift  average  for  any  ensemble  of  dividends 
which  satisfies  the  above  assumption  is  C;  it  does  not  mean  that  with  divisor 
equal  to  D  ,  the  shift  average  is  C  for  each  individual  dividend. 

If  the  distribution  of  the  shift  count,  S,  were  known,  <S>  could  be 
evaluated  by  use  of  the  relation  which  defines  expectation: 


<S>  =  E  m  •  Pr[S  =  m]  (l.l) 

m=l 
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Pr[S  =  m]  can  be  determined  by  converting  f(x),  the  distribution  of  R.  to 
distribution  of  normalized  partial  remainders,  g(x),  using  equation  2:(l.7)j 
then: 


Pr[S  =  m]  =  g(D  +  --)  -  g(D  +  — 7=-)  +  g(D  -  — ~)  -  g(D  -  — ) 
B  v     m         pm+l '       B  v    Qm+1    to  v     m ' 

(1-2) 


Alternatively,  <S>  can  be  determined  without  knowledge  of  the  dis- 
tribution of  S  if  we  make  use  of  a  rather  basic  property  of  ergodic,  aperiodic 
Markov  processes: 


if   M.  =  mean  recurrence  time  of  state  j, 

J 

and  n     =  absolute  probability  of  state  j,  (l-3) 

J 

then   M.  =  (n .  )"  . 
J     J 

Since  S  can  be  regarded  as  the  recurrence  time  of  the  set  of  states 
for  which  x  >  —  (where  steps  of  the  division  process  correspond  to  units  of 
time  in  the  terminology  of  a  Markov  process),  equation  (1.3)  implies: 


<S>"1  =  Pr[R.  >  |] 


or 


<S>  -  - 


1  - 1<£) 


Therefore,  to  find  the  shift  average,  it  is  only  necessary  to  find 
the  area  under  the  density  function,  F(x),  for  x  >  —  (or  x  <  -,  whichever  is 


1 

e.g.,  see  Feller,  Chapt.  XV,  Section  5. 
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more  convenient).  This  simplification  in  evaluating  <S>  is  a  by-product  of 
considering  the  distribution  of  R.,  rather  than  that  of  normalized  partial 
remainders . 


3 •  2  Evaluation  of  <S(P>  for  j  <  D  <  1 

F(x),  for  r-  <  D  <  1,  is  given  by  equation  2:(5.l),  and  we  can  easily 
confirm  Freiman's  result  for  this  region  by  use  of  (1.4): 


<£>  =  — lj-   =  g^   I  <  S  <  1  (£.1) 

2D 

For  ■£  <  D  <  p,  F(x)  is  given  by  equation  2:(5.9)>  from  which  we 
obtain: 


f(|)  =  (2D  -  1)|  +  ^(|  -  2D  +  1) 

=  |,  or  (2.2) 

<S>  =  3,  f or  |  <  D  <  | 

17      5 
The  next  adjacent  interval  in  the  range  of  D  is  — o  ^  D  <  tt,  which  appears  in 

Table  1  as  the  n  =  2  entry  for  sequence  N(n)S.   Since  the  N(2)S  sequence 

involves  three  steps ;,  there  are  four  discontinuity  pairs  in  F(x):  2D  -  1, 

4D  -  2,  8D  -  4,  and  Ikh  -   8,  together  with  the  respective  2D  complements.  The 

area  under  F(x),  x  <  — }    can  be  determined  by  use  of  equation  2:  (5.6).   First: 


F(|)  =  F(D)  +  aT 


+  -L(.  1   )  =  i  ( 


(2.3) 


2D       2DV2T   _      '        2DV2T  _   1' 
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where,  as  defined  in  Chapter  2,  T  is  the  number  of  discontinuity  pairs,  and 
aT  is  the  jump  in  F(x)  at  one  of  the  two  discontinuities  in  the  (p,  2D  -  — ) 
interval,  provided  D  is  not  in  the  set  {D ).   The  area  in  question  is  then 
expressible  as : 

f(|)  =  |  •  P(|)  +  (2D  -  1)&1   +  (4D  -  2)a2  +  (k   -  6D)a3 


[1  +  BD.1  +  iS^  +  i4iiB]        {2.k) 


D(2T  -  1) 


_  2 
3 

(Note  that  2D  -  1,  hL  -   2,  and  k   -  6D  are  the  only  discontinuity  arguments  less 
than  -  for  D  in  this  interval.)  Therefore  from  (1.4): 

<S>  =  3,        for  U  <  D  <|  (2.5) 

17 
Thus  far,  the  shift  averages  for  all  divisors  in  the  range  -~k  <   D  <  1 

have  been  found.   In  addition,  it  can  be  readily  ascertained  from  Figure  h, 

using  equation  (1.4),  that  the  shift  average  for  the  point  D  =  .6  is  also  3. 

17 
The  conjecture  that  <S>  =  3  holds  for  the  interval  .6  <  D  <  ^  is  a  natural 

consequence  of  these  observations.  However,  the  number  of  solution  intervals 
in  this  small  interval  of  D  is  infinite,  with  the  number  of  discontinuities  in 
F(x)  becoming  arbitrarily  large  in  an  infinite  number  of  distinct  neighborhoods, 

The  complex  behavior  of  F(x)  for  ^  <  D  <  4   notwithstanding,  the  con- 
jecture raised  regarding  <S>  in  this  region  can  be  proved  by  use  of  the 
symmetry  property,  derived  in  Section  2.3,  as  demonstrated  by  the  following 
two  theorems . 
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Theorem  3«2a:   There  are  no  discontinuities  in  F(x)  in  the  interval 

h   -  6D  <  x  <  4D  -  2,  for  |  <  D  <  |. 

5      o 


N       N     N  N     N        N 

I      1     I    2   I   3  I  I   3     2   I      1 

2D  -  1       1-D  to-2   J        2D-J       3D  -  1 

2  2 


2  -  2D 


FIGURE  7 


Proof : 

1.  Let  N, ,  Np,  and  N^  denote  the  numbers  of  discontinuity  arguments  appearing 
in  the  intervals  2D-l<x<l-D,  l-D<x<^D-2,  and  ^D  -  2  <  x  <  ^ 
respectively.   Due  to  symmetry,  the  same  numbers  of  discontinuities, 
respectively,  appear  in  3D  -  1  <  x  <  1,  2  -  2D  <  x  <  3D  -  1,  and 
2D--<x<2-2D.   See  Figure  "J.      (Note  that  the  ordering  of  the  end 
points  of  these  intervals  with  respect  to  size  remains  as  shown  in  Figure  7 
for  the  interval  of  D  in  question. ) 

2.  For  D  <  f,    and  2  -  2D  <  x  <  1,  equation  2:(l.5)  reduces  to: 

F(x)=|F(|)  (2.6) 

In  words,  F(x)  for  2  -  2D  <  x  <  1  is  a  copy  of  F(x)  for  1  -  D  <  x  <  .5 
(though  with  horizontal  scale  doubled  and  vertical  scale  halved);  therefore: 


N2  +  N3  -  N2  +  Hx 


or 


(2.7) 


N3  =  N1 
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3.  Due  to  equation  (2.6),  all  N_  discontinuities  in  kD   -  2  <  x  <  —  map  into 
8D  -  h  <  x  <  1.   This,  together  with  equation  (2.7),  implies  no  discontinu- 
ities inl-D<x<4D-2  map  into  3D  -  1  <  x  <  1. 

k.     Therefore,  the  Np  discontinuities  in  1  -  D  <  x  <  ^D  -  2  are  related  to 
those  in  2  -  2D  <  x  <  3D  -  1  in  two  respects:   symmetry,  and  mapping  via 
equation  (2.6).   Due  to  the  latter,  if  x  denotes  a  discontinuity  argument 
such  that  1  -  D  <  x  <  4D  -  2,  the  jump  at  2x  is  half  as  large.  Attempting 
to  select  the  smallest  (nonzero)  jump  among  the  N_  in  1  -  D  <  x  <  kl)  -   2 
would  therefore  lead  to  a  contradiction,  since  another  jump,  half  as  large, 
would  exist  in  2  -  2D  <  x  <  3D  -  1,  and  also  in  the  first  interval  due  to 
symmetry.   Hence,  N  =  0. 

5.   In  step  3,  the  upper  limit  on  one  of  these  symmetric  discontinuity-free 
intervals  was  extended  from  3D  -  1  to  8D  -  h.      The  symmetric  mate  to 
(2  -  2D,  8D  -  k)   is  {k  -   6D,  UD  -  2).   Q.E.D. 

Theorem  3-2b:  <S>  -  3  for  |  <  D  <  |. 
5       8 

Proof ; 

(The  proof  consists  of  relating  areas  under  F(x)  for  certain  intervals 
of  x.   Recognizing  expressions  of  the  form:   f(z  )  -  f(z  ),  as  the  area  under 
F(x),  for  z  <  x  <  z  may  be  found  helpful.) 

1.  Using  equation  2:  (1.4)  for  D  <  ■?:  <, 

f(^D  -  2)  -  f(2D  -  l)  =  f(3D  -  l)  -  f(2D  -  |)  +  f(2D  -  l)  -  f(D  -  |) 

(2.8) 

2.  Since  F(x)  =  |,  for  0  <  x  <  2D  -  1, 

f(2D  -  1)  -  f(D  -  |)  =  1  -  ^  (2.9a) 


and 
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f(2D  -  1)  -  2  -  |  (2.9b) 

3-   Substituting  (2.9)  into  (2.8): 

f(4D  -  2)  -  f(3D  -  1)  -  f(2D  -  |)  +  3  -  ^  (2.10) 

U.   Due  to  symmetry,  and  the  fact  that  F(D)  =  — : 

F(D  +  C)  +  F(D  -  C)  -  i,    for  any  C  <  D  (2.1l) 

Therefore,  the  sum  of  areas  under  F(x)  for  two  symmetrically  located 

W 
intervals,  each  of  width  W,  is  — ,  i.e., 

/Z  Z   ""  z 

2  {F(D  +  y)  +  F(D  -  y)}dy  =  2  p  X 

Zl 

provided  all  arguments  are  in  the  (0,  2D)  interval.   An  application  of  this 
relation  is : 


]  l    D  "  I       1 

f  (±)  -  f(l  -  D)  +  f (3D  -  1)  -  f(2D  -  i)  =  — ^   -I-535     (2*12) 


Substituting  (2.12)  into  (2.10): 

fj(to  -  2)  =  4  -  §  -  [f(|)  -  f(i  -  d)] 


or 


(2.13) 


*(§)  ->  -|-  [f(^D  -  2)  -  f(l  -  D)] 
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5.  From  Theorem  3«2a: 

F(4D  -  2-)  =  F(l  -  D)  (2.l4) 

Using  2: (1.5)  for  D  <  i   and  x  =  4D  -  2- : 


F(4D  -  2-)  =  |[F(2D  -  1-)  +  F(3D  -  1-)] 


=  ^  +  |F(3D-  1-) 


(2.15) 


From  equation  (2.11):: 


F(3D  -  1)  +  F(l  -  D)  =  i  (2.16) 


After  substitutions,  the  last  three  equations  yield: 


F(4D  -  2-)  =  ~  (2.17) 


6.  The  area  under  F(x),  l-D<x<4D-2,  can  now  be  evaluated,  using 
(2.14)  and  (2.17): 


f(kV  -   2)  -  f(l  -  D)  =  ^[k-D  -   2  -  (1  -  D)]  =  f  -   |       (2.18) 


Using  this  result  in  (2.13),  we  obtain: 


f(|)  =  |  (2.19) 


which  produces  the  desired  relation  in  (1.4).  Q.E.D. 
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The  fact  that  the  shift  average  is  equal  to  3  for  all  divisors 

3       3 

in  —  <  D  <  r-  was  indicated  by  use  of  simulation  studies  by  Freiman,  though 

17      3 
proved  only  for  the  subinterval  ^h  <  D  <  t-.   This  result  has  been  used  by 

Professor  G.  Metze  in  proposing  a  modification  to  S-R-T  division.  Briefly,  the 
Metze  form  of  S-R-T  division  uses  a  variable  threshold  value,  rather  than  the 
constant  -^}    to  determine  whether  a  given  step  involves  an  addition  (or  subtrac- 
tion) of  the  divisor.   The  threshold  is  selected  at  the  start  of  any  particular 
division,  based  on  an  inspection  of  the  leading  digits  of  the  divisor;  if  K 
denotes  the  threshold,  then  the  relation  to  be  satisfied  in  order  to  maximize 
the  shift  average  is : 

|  K  <  D  <  |  K 


or 


2  5 
|D<K<7D 

3  6 


3.3  <S(D)>  f or  \  <   D  <  | 

For  D  >  .6,  <S(D)>  has  been  found  to  be  a  relatively  straightforward 
function.   However,  in  the  interval  —  <  D  <  .6,  the  complex  behavior  of  the 
partial  remainder  density  is  reflected  in  <S(D)>  as  well. 

To  evaluate  <S(D)>  for  the  families  of  solution  intervals  for  which 
F(x)  was  determined  in  Chapter  2,  equation  2: (5.6)  will  be  used  to  find  f(^). 
Due  to  symmetry,  together  with  Theorem  2.5b,  one  of  each  pair  of  discontinuity 
arguments  in  less  than  — ,  with  the  exception  of  the  one  pair  in  (— ,  2D  -  — ). 

If  x*  denotes  which  of  the  pair  of  symmetric  discontinuity  arguments,  x^  x^t 
is  less  than  — ,  then: 
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f(|)  ■  \  F(|)  +  Z   a.x*  (3.1) 

i=l 

where  F(p)  is  given  by  equation  (2.3).   As  an  example,  consider  the  infinite 
set  of  intervals  in  the  range  of  D  for  which  the  sequence  of  operations  in 
2:  (5.2)  was  N(n),  n  =  1,  2,  . ..;  viz.: 


n+1  n+1  _ 

5-=  <  D  <  - — = = 

2n+2        2n+2  _  k 


Using  (3«l)  and  2:  (5.6),  and  observing  that  T  =  n  +  1  in  this  case,  we  find 
that: 


1  Pn_1 

f(o)  "  £\ [1  +  n(2D  -  1)],   n  =  1,  2,  ...     (3.2) 

2  D(2n+1  -  1) 


for  this  family  of  intervals.   As  a  second  example,  consider  the  family  of 
intervals  corresponding  to  the  operations  N(n)S  in  2:  (5-2);  T  =  n  +  2,  and  the 
expression  is : 

'<!>  -    n+f     [1  +  "<2D  "  ^  +  2D  ~  ^^  "  ^      (3.3a) 
2    D(2n+2  -  1)  2n 


=  —TZr [1  +  (n  -  1)(2D  -  1)  +  -2-]         (3.3b) 

D(2n+2  -  1)  2n  1 

In  (3.3a),  the  term  n(2D  -  l)  results  from  the  n  multiples  of  2D  -  1  (i.e., 
2J(2D  -  l),  j  =  0,  1,  ...,  n  -  l)  which  are  less  than  -r.      Since  the  sequence 
terminates  after  one  subtraction,  the  discontinuity  argument  2  (2D  -  l)  must 
exceed  2D  -  — ,  and  therefore  its  2D  complement  must  be  less  than  — ,  and  is 
represented  by  the  last  term  in  (3. 3a).   Table  2  is  a  list  of  the  expressions 
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TABLE  2.      EXPRESSIONS  FOR  f(|) 


2:  (5.2)   Sequence 

f(|) 

1)       N(n) 

n  =  2,   3,    ... 

f(|)  "  —|li [1  +  n(2D  -  1)] 

2          D(2n+1  -   1) 

2)       N(n)S 

n  =  2,    3,    ... 

*&  =  ~ nC El  +  (n  -   D(2D  -   1)  +  -2r] 

2          D(2n+2   -   1)                                                     2n  a 

3a)     N(n)SN 

n  -  2,   3,    ... 

1  Pn+1                                                    T) 

f  2    =  -7w5 [1  +  n(2D  "  1}  +  \] 

2  D(2n+3   -   1)                                          2n 

3b)     N(n)SS 

n  -  3,  ^    ... 

t(&  =  — C [1  +■  (n  -  2)(2D  -  1)  +  --?-] 

2          D(2n+3   -   1)                                                       2n  2 

4a)     N(n)SNN 
n  *  3,  fc, 


*(i)  ■ 


,n+2 


>2'       ^n+"5 "  [1  +  (n  +  1)(2D  "  1)] 


D(2ur"   -   1) 


4b)     N(n)SNS 
n  -  2,   3, 


n+2 


f  ri)   =  2 


—   [1  +   (n  -  1)(2D  -   l)+—i-«  D] 
D(2X^   -   1)  2n+1 


4c)     N(n)SSN 
n  -  3,  fc, 


'(J)  ■ 


,n+2 


Nn+ 


D(2U"   -   1) 


[1  +   (n  -  1)(2D  -   1)  +  -^J  «D] 


4d)     N(n)SSS 
n  -  *,  5, 


n+2  - 

f(|)  =  ha [I  +  (»  -  3) (2D  -  i)  +  -=s  -D] 

2         D(2n+4  -  1)  2n  X 


-63- 

for  f(p)  for  all  subintervals  of  I  ,  n  =  2,  3;  . ..,  in  which  there  are  n  +  h 
discontinuity  pairs  or  less;  the  interval  of  D  to  which  each  applies  is  given 
by  Table  1. 

In  Table  3>  expressions  derived  for  Table  2  are  used  to  evaluate 
<S(D)>,  using  equation  (1.4),  and  with  n  =  2,  3>  4,  5-   One  interesting 
property  of  <S(D)>  in  this  interval  is  the  appearance  of  "plateaus"  for  certain 
intervals  of  D  (e.g.,  N(2)S,  N(3)SS,  N(4)SSS,  ...).   The  characteristic  which 
identifies  plateau  regions  is  evident  from  examination  of  expressions  for  f (p); 
any  sequence  which  involves  a  total  of  m  normalizations  and  m  subtractions 
yields  an  expression  of  the  form: 

T-2 

f(|)  =  % [1  +  (m.  -  m  )(2D  -  l)  +  Q    •  D]       (3-4) 

2    D(2T  -  1)        1  2 

where  9   is  a  constant.   (Equation  (3«l)  can  be  used  to  confirm  readily  that 

(3.4)  is  true  for  the  general  case: 

a)   If  the  step  required  to  generate  x   ,  given  x  ,  is  a 

subtraction,  then  x  >  2D  -  —  and  the  a   *  x*  product 

K         2  K     K 

in  (3.l)  is  of  the  form 
JC-1 


2D  -  2    (2D  -  1)  +  ... 
2K-1 


2T-2 


D(2T  -  1) 


b)  Analogously,  if  the  kth  step  is  a  normalization,  then 

x  must  be  less  than  — ,  and  the  a   *  x.*  product 
K  2  K     K 

contributes  a  positive  multiple  of  2D  -  1. ) 
Therefore,  when  m  exceeds  m  by  (exactly)  one,  (3-^)  yields  a  constant. 
Conversely,  f (—),  and  therefore  <S(D)>,  can  be  constant  over  an  interval  of  D 
only  when  m  -  m  =  1.   Since  nothing  has  been  assumed  regarding  the  value  of 
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TABLE  3-   EXAMPLE  EXPRESSIONS  FOR  <S>,  USING  TABLE  2  AND  EQUATION  3: (1.4) 


Sequence 

<S> 

in  2:(5.2) 

n  =  2 

n  =  3 

n  =  k 

n  =  5 

l)   N(n) 

7D 
2  -  D 

15D 
8  -  9D 

31D 
2k   -  33D 

6  3D 
6k  -   97D 

2)   N(n)S 

3 

31D 

8  -  3D 

63D 
32  -  35D 

127D 

96  -  131D 

3a)  N(n)SN 

31D 
8  -  3D 

63D 
32  -  35D 

127D 
96  -  131D 

255D 

256  -  387D 

3b)  N(n)SS 

63 
23 

12  7D 
32  -  9D 

255D 

128  -  137D 

4a)  N(n)SNN 

12  7D 

96  -  129D 

255D 

511D 

256  -  385D 

6UO  -  1025D 

4b)  N(n)SNS 

3 

127D 
32  -  IID 

255D 

511D 

128  -  139D 

384  -  523D 

kc)     N(n)SSN 

127D 
32  -  IID 

255D 

511D 

128  -  139D 

384  -  523D 

4d)  N(n)SSS 

255 
103 

511D 

128  -  25D 
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3      3 
D,  all  solution  subintervals  in  the  —  <  D  <  t-  region  must  correspond  to 

sequences  for  which  this  difference  is  unity.   A  second  consequence  is  that 

F(x)  must  have  an  even  number  of  discontinuity  pairs  for  intervals  of  D  in  which 

<S(D)>  is  constant,  since  the  number  of  operations  in  2: (5 .2)  is  odd  if 

mi  "  m2  =  lo 

In  principle,  Table  3  could  be  extended  indefinitely  to  obtain 

algebraic  expressions  for  <S(d)>  for  as  great  a  fraction  of  the  „5  <  D  <  .6 

interval  as  we  desire,  although  the  size  of  solution  intervals  decreases  rapidly 

with  larger  numbers  of  discontinuities  and  are  particularly  minute  in  the  lower 

portion  of  this  interval  of  D„   (See  Table  1  for  n  =  3,  ^  and  %)     Alternatively, 

by  evaluating  <S(D)>  numerically,  for  a  sufficiently  fine  grid  of  points  in 

this  interval,  a  graphical  representation  can  be  obtained  to  within  any  desired 

accuracy,, 

Figure  8  was  obtained  with  the  aid  of  ILLIAC  II,  using  increments  of 

.001  for  .5  <  D  <  060   The  discontinuity  arguments,  x*,  were  determined  for 

i  <  30o   By  rewriting  equation  (3»l)  as: 


1  0T-2    r    T-2  r  t  . 

2  D(2T  -  1)  L     i=0  L  2       1+1 

we  observe  that  the  maximum  error  in  f (^)  which  can  occur  as  a  result  of 

-30 
T  <  30,  where,  in  fact,  T  may  be  infinite,  is  of  the  order  of  2 

Using  increments  of  10   and  the  trapezoidal  rule,  a  numerical 

approximation  to  the  integral 

/    <S(D)>  dD 


(3.5) 
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was  found  to  be  2.61917  +  5  x  10   .   This  corresponds  to  the  average  value  of 
<S>  for  D  in  this  range,  if  all  values  of  D  are  assumed  equally  likely.   Using 
this  result,  together  with  3:(lA)  and  theorem  3-2b,  the  average  <S>  for  the 
entire  divisor  range  is  2.6170  for  S-R-T  division. 

Because  of  the  short  computation  time  required  per  point  (less  than 
a  millisecond),  it  was  found  practical  to  investigate  in  detail  the  behavior 
of  <S(D)>  for  certain  portions  of  the  divisor  range.   The  following  are  listed 
as  empirical  results : 

a)  The  left  limit  of  the  plateau  region  in  <S>  which  has 

9  ikk 

D  =  -rz   as  the  right  end  point  is  D  =  37^. 

b)  The  above  plateau  is  a  local  minimum  in  <S>.   Moving 

Ikk 
to  the  left  of  D  =  7^j}   <S>  has  numerous  local  maxima, 

each  of  which  is  a  plateau  region.   The  limiting 

2k  96 

maxima  is  a  point,  viz.,  D  =  r-r-,  at  which  <S>  =  ^=-. 

While  the  difference  of  (approximately)  0.01  in  the 

2k  ikk  9 

values  of  <S>  at  D  =  r-r  compared  to  the  -^p=  <   D  <  y^ 

plateau  may  have  little  practical  significance,  this 

behavior  is  sufficiently  contrary  to  intuition  to 

warrant  future  investigation.   Noteworthy  is  the  much 

larger  peak  which  Freiman's  simulation  of  the  so-called 

9 
. 75-1 • 0-1. 5  method  yields  for  D  near  y£. 

c)  The  nonmonotomic  behavior  of  <S>  described  in  (b) 

17 
occurs  again  near  D  =  — ,  and  a  plausible  conjecture 

is  that  it  will  be  found  near  all  D  =  —  +  (p)  , 

n  =  k,    5>  ...,    due  to  the  similarities  in  these 

neighborhoods  noted  in  Chapter  2. 


h.      SUMMARY  AND  CONCLUSIONS 

4.1  Summary 

Treating  digital  division  as  a  Markov  process  in  order  to  provide  a 
basis  for  evaluating  and  comparing  division  algorithms  was  the  contribution  of 
Dr.  C.  V.  Freiman.   If  the  set  of  dividends  encountered  is  assumed  to  have  some 
degree  of  randomness,  so  that  it  may  be  approximated  by  a  piecewise  uniform 
probability  distribution,  then,  as  demonstrated  in  Chapter  1,  the  distribution 
function  describing  partial  remainder  magnitudes  approaches  a  limit  which  is 
determined  only  by  the  divisor,  D,  and  in  particular,  is  independent  of  both  the 
dividend  distribution  and  the  step  numbers  in  the  decision  process.   This  stationary 
distribution  is  itself  piecewise  uniform,  and  the  intervals  of  uniformity  are 
the  states  of  the  Markov  chain.   However,  a  straightforward  application  of  the 
theory  of  finite  Markov  chains  in  the  case  of  S-R-T  division  was  found  imprac- 
ticable, for  the  following  reasons : 

a)  For  S-R-T  division,  the  number  of  discontinuities  in 
the  partial  remainder  stationary  density  function 
becomes  indefinitely  large  in  an  infinite  number  of 
distinct  neighborhoods  in  the  divisor  range,  as 
demonstrated  in  Chapter  2. 

b)  When  solving  for  the  probabilities  of  the  individual 
states,  or  intervals,  the  number  of  unknowns  in  the 
linear  equation  set  must  at  least  equal,  and,  in 
general,  will  exceed  the  number  of  discontinuities 
in  the  density  function. 

c )  Solving  one  set  of  linear  equations  yields  the  partial 
remainder  distribution  for  only  one  value  of  the 
divisor. 
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As  an  alternative  to  the  finite  Markov  chain  approach,  a  functional 
relation  was  derived  which  was  shown  to  be  a  necessary  and  sufficient  condition 
for  the  stationary  density  function  (2:(l.5)).   By  considering  the  Fourier 
series  representation  of  a  periodic  extension  of  the  density  function,  and 
applying  the  functional  relation  just  mentioned  together  with  the  knowledge  that 
the  total  variation  (per  period)  is  bounded,  it  is  proved  that  the  density 
functions  for  all  values  of  D  must  satisfy  an  important  symmetry  relation 
(Section  2.3).   This  symmetry  is  used  to  derive  a  number  of  other  properties 
which  lead  to  a  straightforward  algorithm  (2:  (5.2))  for  determining  all  discon- 
tinuity arguments  in  the  density  function,  and  an  equation  for  the  size  of  the 
respective  jumps  (2: (5.6)).   The  limiting  distributions  being  sought  may  be 
expressed  in  terms  of  these  jumps.   The  solutions  thereby  obtained  apply  for 
intervals  in  the  range  of  D,  and  by  grouping  these  intervals  on  the  basis  of 
the  sequence  of  steps  used  in  2:  (5.2),  infinite  families  of  intervals  are 
formed.   Properties  of  the  array  of  solution  intervals  are  derived  by  consid- 
ering, among  other  things,  the  relations  which  must  exist  between  adjacent 
intervals . 

While  the  partial  remainder  stationary  density  function  is  interesting 
in  itself,  a  practical  application  lies  in  evaluating  the  shift  average,  <S>, 
which  is  the  average  number  of  quotient  bits  generated  per  iteration;  this 
problem  is  considered  in  Chapter  3.   Since  the  shift  average  is  equal  to  the 
mean  recurrence  time  of  partial  remainders  which  exceed  0.5  in  magnitude,  <S> 
is  found  to  be  simply  related  to  the  partial  remainder  distribution  function 
(3:(l.*0).   Properties  derived  in  Chapter  2  enable  proving  that  <S>  =  3  for  the 
entire  interval,  0.6  <  D  <  0.75^  a  result  which  was  previously  not  obtainable 
due  to  the  complex  behavior  of  the  density  function  for  D  near  0.6.   For 
0.5  <  D  <  0.6,  the  complicated  behavior  of  the  partial  remainder  density 
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function  is  reflected  in  <S>  (Table  3),  and  relations  derived  in  Chapter  3  are 
used  to  evaluate  <S>  with  aid  of  ILLIAC  II,  for  D  in  that  interval  (see 
Figure  8). 

These  results,  together  with  equation  3:(2.l),  the  expression  for 

3 
<S>  when  r-  <  D  <  1,  satisfy  a  goal  defined  at  the  outset:   determination  of 

the  shift  average  for  the  entire  range  of  divisors. 

h.2      Conclusions  and  Areas  for  Future  Investigation 

In  recent  years  there  have  been  a  number  of  efforts  to  better  under- 
stand the  digital  division  process  when  more  than  a  minimum  number  of  divisor 
multiples  are  available,  and  to  apply  this  insight  in  formulating  algorithms 
which  reduce  the  number  of  nonzero  quotient  digits .   Usually  the  avenue  of 
approach  is  to  consider  bit  patterns  which  can  occur,  then  to  select  a  set  of 
rules  which  provides  best  performance  in  each  example.   This  is  often  an 
iterative  process. 

Knowledge  of  stationary  distributions  of  partial  remainders  for  the 
various  classes  of  division  algorithms  provides  a  valuable  theoretical 
foundation  in  this  search  for  optimum  algorithms.   Thus,  for  example,  the  shift 
average  obtained  for  S-R-T  division  is  applicable  to  an  entire  class  of  what 
might  be  called  single-threshold  division  processes,  which  can  be  described  as 
follows : 

a)  the  divisor  multiples  available  are:  mD,  0,  and  -mD, 
where  m  >   1, 

b)  there  is  a  positive  number,  K,  such  that  K  <  |mD|  <  2K 
is  true,  and 

c)  K,  which  in  the  general  case  is  a  function  of  D,  is 

used  as  the  threshold  level  with  which  each  R.  is 

i 
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compared  to  decide  whether  or  not  the  divisor  is  to 

be  used  in  forming  R. 

&  i+l 

The  abscissas  in  the  <S(D)>  versus  D  relation  derived  in  Chapter  3  may  be 
regarded  as  —  D  rather  than  D.   Metze's  modification  of  S-R-T  division  uses 

CIS. 

this  change  of  scale  effected  by  varying  K  to  take  advantage  of  the  maximum 
which  <S>  was  found  to  assume  when  0.6  <  D  <  0.75  is  true. 

An  important  question  then  is:   can  any  similar  generalizations  be 
made  to  optimize  multiple -threshold  algorithm,  perhaps  using  a  property 
analogous  to  the  relative  independence  of  m  and  K  in  the  single-threshold  case? 
(Only  the  inequalities  in  (b )  above  must  be  satisfied.)   In  other  words,  when 
more  than  one  (nonzero)  multiple  of  the  divisor  magnitude  is  assumed  to  exist, 
to  what  extent  can  the  solution  for  a  specific  algorithm  be  generalized  in 
order  to  aid  in  formulating  an  optimum  one? 

Derivation  of  partial  remainder  distributions  for  any  useful  two- 
threshold  algorithms  (MacSorley1 s  paper  includes  examples  of  such  algorithms) 
is  complicated  by  the  fact  that  a  number  of  ranges  of  the  divisor  must  be 
analyzed  individually,  since  the  thresholds  vary  with  D.   However,  the  observa- 
tions that  a)  the  piecewise  constant  functions  which  describe  the  stationary 
density  can  be  formed  by  locating  discontinuity  arguments  as  an  ordered  process, 
starting  with  discontinuities  in  the  defining  equation,  and  b)  the  size  of 
jumps  at  successive  arguments  in  this  process  form  a  geometric  progression, 
should  expedite  these  derivations. 

In  addition  to  considering  binary  division  with  further  redundancy 
added,  there  is,  of  course,  the  problem  of  considering  higher  radices. 

A  second  area  where  further  work  is  needed  is  the  modification  of  the 
<S(D)>  versus  D  relation  to  account  for  the  fact  that,  in  parallel  computers, 
shifting  R.  an  arbitrary  number  of  places  in  one  step  is  not  possible.  One 
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approach  to  this  problem  would  be  that  outlined  by  equations  3:(l.l),  (1.2), 
and  (1.3)?  with  the  index,  m,  in  (l.l)  confined  to  a  finite  range  corresponding 
to  the  shift  paths  available. 

Finally,  a  more  thorough  understanding  of  some  peculiarities  in  the 
stationary  distribution  of  R.  for  S-R-T  division  would  be  interesting  if  not 
useful.  A  few  of  the  questions  which  remain  are: 

a)  Can  the  nonmonotonic  behavior  of  <S(D)>  for  0.5  <  D  <  0,6 
be  related  to  divisor  bit  configurations?  The  fact  that 
"bad  examples,"  i.e.,  cases  for  which  S-R-T  type  division 
is  clearly  nonoptimal,  often  involve  divisors  of  the 

1       In 

form  —  +  (— )  ,  and  that  the  local  minima  in  <S(D)>, 
which  are  plateau  regions,  include  these  points,  indicates 
an  intuitive  explanation  for  these  minima  and  what 
determined  their  end  points  may  not  be  difficult. 

b)  What  significance  can  be  attached  to  the  accumulation 
points  of  solution  intervals,  i.e.,  neighborhoods  of 
D  which  lead  to  distributions  of  R.  with  arbitrarily 
large  numbers  of  discontinuities? 
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